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/  Abstract 

^  ;  “-rr- - 

Application  of'lolal  \  ariation  Diminishing  (T\  D")  schemes  to  both  inviscicl 
and  viscous  flows  is  considered.  The  mathematical  and  physical  basjs  of  TVD 
schemes ^s .discussed.  First  and  second-order  accurate  TVD  schemes/and  a  second- 
order  accurate  Lax-W'cndrofr  scheme,,  are  tised  to  compute  solutions  to  the  Riemanii 
problem  in  order  to  investigate  the  capability  of  each  to  resolve  shocks,  rarefactions, 
and  contact  surfaces.  .Second  order  finite-volume  and  finite-difference  TVD  schemes 


are  used  to  obtain  solutions  to  iiniscid  supcr.sonic  and  transonic  cascade  flow  piob- 
lems.  TVD  schemes  arc  shown  to  be  superior  to  the  Lax-VVendroff  family  of  scheme's 
for  both  transient  and  steady-state  compulations. 


T\'D  methodology  is  extended  to  Vhe  solution  of  viscous  flow  problems.  .A  first- 
order  time  accurate,  second-oruei  space  accurate  algorithm  is  contrasted  against 
a  second-order  time  and  space  accurate  algorithm  for  khe  solution  of  the  viscous 

Burgers’  equation.  .\'ecc.ssity  of  using  the  fully  second-order  accurate  algorithm  at 

•  • 

low  Reynolds  numbcr.s  is  shown.  .Solutions  aie  computed  tolhe  problems  of  laminai 
shock-boundary- layer  interaction  and  unsteady,  laminar,  shock-induced  heat  tiansfei 
insing  (he  new  algorithms.  These  algorithms  provide  the  capabilityjjor  the  first  time. 
'  to  accurately  predict  separation,  reattac.hment.  and  pressure  and  akin  friction  piofdes 
for  shock  boundary-layer  inleraction.  .AdditionalhVextiemely  accurate  coinpaiisoii 
with  theory  and  experiment  is^evident  for  the  unsteady,  shock-induced  heal  tiansfei 
problem. /riiese  solutions  are  lont lasted  against  solutions  computed  with  the  Be.  m- 
VVarming/algorithm.  and  the  T\'D  solutions  are  shown  to  be  vastly  superior. 


HIGH-RESOLUTION  TVD  SCHEMES  FOR  THE  ANALYSIS  OF 
I.  INVISCID  SUPERSONIC  AND  TRANSONIC  FLOWS 
II.  VISCOUS  FLOWS  WITH  SHOCK-INDUCED 
SEPARATION  AND  HE.AT  'niANSFER 


I.  Introduction  to  Part  I 


LI  Overview  of  Pari  I 

Part  I  begins  with  a  historical  look -at  the  (levelopr''ent  of  what  has  become 
known  as  the  Total  Variation  Dimini.'jhing  (T\’D)  class  ol  schemes  lor  .-'olving  hypej- 
bolic  conservation  laws.  Conditions  nece.ssai\  fur  a  finite-PilTcrenre  scheme  to  yield 
physically  meaningful  solutions  arc  distus.sed.  l)e\elopmei*t  of  a  .second-order  accu¬ 
rate  TVD  scheme  for  scalar  conservation  laws  is  detailed,  along  with  the  means  of 
extending  il  to  .systems  such  as  the  Eulei  equations  of  ga.sdynamics.  .A  brief  discu.s- 
sion  of  the  Euler  equations  is  undertaken  in  the  context  of  applying  T\'D  schemes  to 
their  solution.  First-order  T\’D.  second-order  1'VD.  and  .second-oider  Lax-VVendrolf 
schemes  are  applied  to  the  Riemann  prol'lent  ti»  determine  ilie  al)ility  of  eai  h  to 
re.solve  the  relevant  features.  Two  second-oider  T\’D  schemes  foi  .solving  .systems  of 
equations  in  two  space  dimensions  <tre  coveied.  'Phese  two  schemes  aie  then  applied 
to  the  .solution  of  both  supersonic  and  tran.soni*  (a.--  idr  flow  problems.  J’VD  algo¬ 
rithms  are  shown  to  be  vastly  superior  to  the  Lax Wendiolf  fc  ;  of  algorithms  for 
both  transient  and  steady-state  solutions. 


1.3  The  Gtnc$i$  of  TVD 

Total  Variation  Diminishing  (T\'D)  schemes,  originally  refered  to  as  'Jbtal  Vai  i- 
ation  Nonincreasing  (T\'’NI).  first  appeared  m  1!)S3  with  the  publication  of  Marten  s 
[figli  Resolulion  Schemas  for  Ifgpcrholit  Con.-ifrralion  /,</ic.s  [22).  In  general.  T\’f) 
.scheme.^  are  arrived  at  by  applying  a  first  order  accurate  numerical  method  to  an 


'ippropriately  modified  flux  function  thus  yielding  a  method  vnat  is  second-oidci  ac¬ 
curate  except  near  points  of  extrema  of  tlie  solution.  The  genesis  of  the  T\’D  cla.ss 
ol  finite-difference  schemes  can  be  traced  to  1976  \vh(  n  llarten,  Hyman,  and  Lax  au¬ 
thored  On  Finitc-DiJJerence  Approxirndiions  and  Eni  o  *  ''nndilioni.  for  Sliock-i.[2'^]. 
This  work  first  addressed  the  question  of  whethc,  xn  •-.fifference  approximations 
to  the  solution  of  hyperbolic  conservation  laws  ^  .. .e;-.  'c  i-he  ph.ysically  relevant 
solution.  This  is  of  interest  because  w'eak  solutions  to  -  ch  conservation  laws  are  not 
uniquely  determined  by  initial  values,  but  require  an  cnitiODy  condition  be  met  to 
converge  to  the  particular  physical  solution  (2.3). 

In  the  mid  i970’s,  llarten  w'as  also  working  on  his  artificial  co;npie.s.sion  method 
(.‘VCM)  [19]  to  modify  standard  finite-difference  schemes  in  an  effort  to  prevent  the 
smearing  of  contact  surfaces  and  improve  shock  resolution  (20.  21 ).  Prioi  lo  this 
effort  Hartcn  states  that  the  standard  finite-difference  schemes  in  use  tvpicallv 
smeared  shocks  over  3-5  cells  while  the  width  of  the  contact  surface  behaved  as 
jg  {^otal  number  of  time  steps  taken  and  R  is  the  sthenic's  order 
of  accuracy.  Marten’s  .4CM  also  addressed  the  fact  tliat  schemes  of  orclci  greate: 
than  one  produced  overshoots  and  undershoots  around  the  discontinuity  [20j  i.nd 
forced  the  approximated  solution  lo  be  nonphysical  [23].  Marten's  .\CM  modifica¬ 
tions  to  existing  schemes  provided  the  1  mndation  foi  the  new  class  of  TV’D  schemes 
prevsented  in  his  1983  paper. 

The  rigorous  mathematical  foundation  of  'Pn’D  schemes  is  mainly  confined  lo 
scalar  linear  and  nonlinear  conservation  laws  and  is  painstakingly  oullineil  in  r«>l- 
erences  [23]  and  [22].  Computational  fluid  dynamicists  are  interested  in  applying 
TVD  .schemes  lo  .systems  of  nonlinear  hyperbolic  conservation  k.w.s.  such  a.*-  the  l‘lu- 
ler  tqualions  of  gasdynamics.  Therefore.  Marten  details  the  application  of  TVD 
methodology  to  I-D  systems  using  Roe's  appioximate  Riemann  solvei  <md  piovido 
an  example  of  its  extension  lo  2-D  using  .Strang's  dimensional  splitt  ing  [22].  The  oi  ig- 
inai  Marten  scheme  u'as  a  second-ordei  acriirale  explicit  method  but  wa.s  c.\tended 
lo  a  .second-order  accurate  implicit  method  by  ’^’ce  and  Marten  [-13]. 

The  high-resolution  TV’D  approach  soon  gathered  favor;  explicit  and  implicit 
variations  were  then  applied  to  the  Ruler  equations  in  general  gcomcliies  by  Y<’e  tiiul 
Kutler  [14]  and  by  Yee  and  Marten  [16].  Later.  Wang  and  Widiiopf  furthei  extended 
Marten's  TVD  me.thodology  to  a.  finite-volume  S(  heim  for  the  Ruler  equations  [l8j. 


TVD  algorithms  have  continued  to  develop  over  the  past  decade,  llarten's  origi¬ 
nal  scheme  was  of  the  upwind  variety,  meaning  that  the  modification.''  to  the  flu.\ 
function  are  applied  base’l  on  the  direction  of  wave  propagation,  or  chara*  i-'rist  ic 
d'rection.  Symmetric  algorithms  have  since  come  into  use  wliore  the  modilif  ation-' 
are  applied  without  rega'-d  to  the  characteristic  directions.  MetluKls  are  also  avail¬ 
able  for  partial  differential  equations  with  source  terms  and  stiif  source  term.s.  ^'e(‘^s 
1989  publication, .4  Class  of  Uigh-Ile^olulion  Exr'icii  and  Implicit  Shock-Caphivnuj 
Methods  [45],  piv  ■  ides  derailed  information  on  numerous  versions  of  d'A'D  algorithm.s 
and  e.xamples  of  their  application  to  numerous  problems. 

I.S  Hyperbolic  Conservation  Laws  and  'j  \  D  Methodology 

The  present  section  provides  a  description  of  the  hyperbolic  cou.''ervation  law.s 
for  which  TVD  schenies  provide  solutions.  The  rec|U'rements  for  uniqueness  of  a 
solution  to  the  initial  value  problem  are  given  along  with  the  necessary  conditions 
to  guarantee  convergence  of  a  finite  difference  approximation  to  this  solution.  A 
summary  is  provided  of  the  methodology  behind  the  construction  of  Marten’s  original 
second-order  accurate  TVD  scheme. 

1.3. 1  Finile-DiJJcrence  Schemes  ana  Oleinik's  Entropy  Condition.  The  ;)t<  wtil 
analysis  is  concerned  with  weak  .solutions  of  the  initial  value  problem 

«/  -i-  =  0 

—  eo  <  .r  <  CO  (1.1) 

t/(.r,0)  =  o{.r) 

where  u(x,t)  is  a  column  vector  of  m  unknowns,  /(«)  i,.;  the  flux  vector  of  in 
components,  and  e>(.r)  is  the  initial  data.  Mci  1.1  is  hyperbolic  if  all  eigenvalue'' 
«*(u), - (t”’(u)  of  Mie  .Jac.obian  matrix 

■4(  «)  =  ./„  (l.’i) 

are  real  and  the  set  of  right  eigenvector.''  /?'(») . /?’"(»)  is  complete  [22]  over  ilu' 

domain. 


Following  Marten  [22].  con.'>i(ler  sy.st.ein.',  o!  conservation  law.s.  IL]  l.l.  po.s.sess- 
ing  an  entropy  rnnclion  l'(u)  defined  such  that 


/  >  0 

I'ufu  =  F„ 


(1.3) 


where  /•’  is  a  I’uncLion  known  as  the  entro|)y  (in.':  [22]. 

The  class  of  all  weak  solutions  to  Eq  l.l  is  too  large  in  that  the  initial  value 
problem  is  not  uniciuc  [23].  .\n  additional  constraining  relation  is  needed  if  the 
scheme  is  to  choo.se  the  phy.‘'ically  relevant  .solution.  This  additional  constraint  is 
known  as  Oleinik's  entropy  condition  and  can  be  e.xpressed  as  [23] 


(.  {u)i  +  /' (</ '.r  S.  d 


(l.-l) 


Let  us  now  consider  numerical  .solutions  to  Eq  1.1  obtained  using  a  (2/.;+  1) 
point  explicit  scheme  in  conservation  I'onn  [23].  A  scheme  is  in  conservation  form  if 
it  can  be  expres.sed  as 


.'i-M 

/ 


(l.o) 


wliere 


fU 


\/i 


(!.()) 


and  A  =  A//A.r.  In  Eqs  l.-o  and  l.fi.  /  is  the  "numerical  "  .  or  mesh,  flux  function 

consisrenl  with  /(«;)  in  that  /(» . u)  =  f{ti)  I'lie  solmion  it  is  apj)roximaled  on 

the  me.sh  by  r"  =  r(jA.r.nA/).  The  nmnericai  .sc  heiiu*  given  by  E(|  l.-a  i.s  consistent 
wiih  the  entropy  condition.  Eq  l.  l.  if 


r"-' 

j 


I  I.’U 


(i.Ti 


where  f"  =  U  (e").  . F  i.s  the  numerical  entropy 

flux  consistent  with  F(j/)sucii  that  F(u . u)  =  F{u) 

The  question  of  convergein e  of  the  finite  clifference  scheme.  Et|  i.-o.  to  the  ap¬ 
propriate  weak  .solution  of  Eq  l.l  mu.st  now  be  addre.ssofl.  The  heme  under  ( on.sid 


I 


eration  is  nonlinear,  so  stability  of  <i  consistent  .scheme  does  not  imply  coineigeiice. 
Hai  ten  [■22]  outlines  three  conditions  uhich.  when  satisfied,  en.sure  coiiveigcnce. 

(1)  The  total  variation  (7’\')  of  (he  finite  difference  scheme  is  uniformly 
bounded,  where 

■/■'■(<■)=  y;  I  (I.S) 

J=-oc 

(2)  The  scheme  is  consistent,  as  A.r  0  .  wuth  Oleinik's  entropy 
condition  for  all  entropy  functions  of  Eq  1.1. 


(3)  Oleinik's  entropy  condition  implies  a  unique  solution  of  the 
initial  value  problem  for  Eq  1.1. 

The  reader  is  referred  to  the  references  given  by  Marten  [22]  for  the  arguincnt.s 
that  imply  convergence  gi\en  satisfaction  of  the  above  criteria.  For  the  present 
work,  the  validity  of  these  criteria  will  be  assumed  and  the  elTort  concentrated  on 
demonstrating  the  development  of  a  .siheme  that  satisfies  criteria  (I)  and  (2)  when 
given  the  third  criterion. 


1.3.2  Development  o(  lInrU  u  s  St  coml-Ordc.r  Senior  TVD  Scheme.  Marten's 
second-ordei  accurate  T\’D  scheme  is  the  product  of  a  nono.scillatory.  firsl-ordci 
accurate  scheme  applied  to  an  appropiiat<'U  modified  llu.x  function  [22j.  Thi.s  .section 
describes  the  properties  of  the  first  ordei  .m  heine  and  outlines  the  piocediire  u.M*d  In 
Marten  to  arrive  at  the  appropriate  modified  fhi.x. 

Consider  the  initial  value  problem  for  a  .scalar  con.servation  law: 

III  -r  /(h)j.  S  II,  4-  iiln}iij.  =  0 

rt(u)  =  /u  -cc  <  .r  <  cc  (I--)) 

ii(x.  0)  =  c>(.r) 


where  o(.r)  is  of  bounded  total  variation.  IMgoious  analy.sis  is  re.slricted  to  the  .valai 
ca.se  bi'cau.se  T\  D  .schemes  are  not  defineti  foi  systems  of  of  non-linear  (on.sercation 


Taws  where  the  spatial  total  variation  of  the  solution  may  increase  clue  to  wave 
interaction  [22). 

A  weak  solution  of  Eq  1.9  has  a  monotonicity  property  (22),  as  a  function  of 
time,  defined  as: 


(1)  No  new  local  e.x'trema  in  .r  may  be  created. 

(2)  A  local  minimum  is  nondecreasing  and  a  local  ma.x'imum 
is  nonincreasing. 


The  monotonicity  property  implies  that  the  total  variation  in  .r  is  nonincreasing  in 
lime.  TV’ (w  (/.'>))  <  7'V  (n  (/| )). 

.\n  explicit.  {21:  +  1),  point  finite- difference  .scheme  in  conservation  form.  a.s 
given  by  Eq  1.5  and  applied  to  Eq  1.9.  can  be  written  as 


"V;-* . 


(I.IO) 


or  in  operator  notation  as 


e"  ■ '  =  L  ■  f” 


(1.11) 


The  scheme  given  by  Eq  I.IO  is  T\’D  if.  for  all  v  of  bounded  total  variation 


TV(  L  •  r)  <  7'V'{e) 


(1-12) 


where  the  total  variation  is  defined  by  Eq  1.8.  Eq  l.ll  represents  a  monotonicity 
prc.serving  .scheme  if  the  operator  L  is  monotonicity  preserving.  That  i.-.  if  r  i."  a 
monotonic  mesh  funct  ion  so  is  L  ■  r.  Tire  hclieme  is  monotone  if  II  is  a  monoit.-nic 
nondecreasing  function  of  each  of  its  21:  -f  1  aigumenis  (2.'lj: 

//„.,(  <r_jt - .'cr:)>0  (l.Iaf 


for  all  i  such  that  —k  <  i  <  k. 


r, 


An  example  of  a  scheme  that  is  not  monotone  is  the  scconcUorcler  accuiate 
Lax-VVendroff  scheme  with 


where  =  fjj-i  —  Vj  .  Therefore  the  discrete  equation  is 


(I.I-!) 


't 


Taking  the  derivative  of  IJ  with  respect  to  tlie  argument  e'T,  yields 


II 


(1.16) 


where  //  =  a\  .  Only  the  case  0  <  //  <  I  need  be  examined  since  the  Law-Wendroff 
scheme  is  unstable  for  u  >  i.  and  Lax-WendrolT  provides  the  exact  .solution  for  //  =  1. 
Clearly,  the  Lax-VVendroff  scheme  is  not  monotone  for  any  ;/  <  I  .  Additionally, 
the  numerical  results  of  reference  {23j  show  tliat  the  Lax-WendrolT  ^chcme  is  not 
monoionicity  preserving. 

The  first-order  accurate  Roe  .«che!ne  provide.s  an  example  of  monotone  behav¬ 
ior.  The  numerical  flux  for  the  Roe  scheme  is 


(i.lTf 


giving  the  discrete  equation 


••f'  =  [/  -  f  -  !"i  {'■:»  -  a-: + 


f  l.lv 


Taking  derivatives  of-  II  u-itii  respect  to  Ccich  of  its  arguments  gives 


IL,.,  = 

=  I-//  {1. 10) 

=  0 

Thus,  II  is  a  moiiolonic,  non-decreasing  function  of  each  of  its  .irguments  >liouing 
that  the  Roe  scheme  is  indeed  monotonic. 

Let  5,\/,  St\  P-,  and  S;wp  denote  monotone.  TV’D.  and  monolonicity  preserving 
.schemes  respectively.  Theorem  2.1  of  reference  (22)  provides  the  hieiarchv  of  the.se 
properties: 

S.\i  C  Stvd  C  S;iii>  (1-'20J 

Thus,  the  Roe  scheme  is  also  TVD  and  monolonicity  preserving. 

A  scheme  in  the  conservation  form  of  Rq  1 .10  that  is  monotone  with  v”  converg¬ 
ing  boundedly  almost  everywhere  to  some  function  u(x.l)  has  two  further  desirable 
properties.  The  theorem  of  Lax*  and  VVendrofF  a.s  given  by  reference  [23]  states  that 
if  the  scheme  is  in  conservation  form  with  v(x.l)  converging  almost  everwhere  to 
M(.f. /).  then  u(x.t)  is  a  weak  .solution  of  Eq  1.9.  The  theorem  of  llarten.  Hyman, 
and  La.\  (23)  stales  that  if  the  scheme  is  jiionotone  in  addition  to  meeting  the  crite¬ 
ria  of  the  Lax  Wendroif  theorem,  then  Oiciiiik’s  entropv  condition  is  .s.itisfied  for  all 
di.sconlinuilies  of  u.  Thus  a  monotone  scheme  .'>.vtisfics  the  convergence  ciilcria  foi  a 
unique  solution  of  the  initial  value  problem  as  staled  in  the  Section  1.3.1. 

.Attention  is  now  focused  on  how  the  properties  of  a  monotone  scheme  are  hel|)- 
fu!  in  constructing  llarten 's  second  order  T\'D  .-ciieme.  Ilailen  .slates  that  monotone 
.schemes  provide  second-order  accurate  .sulu lions  to  the  modified  Eq  (22) 


M}  -f  /(h1-  =  Jv/  '  ->1  u.  Alu,-), 

(1.211 

^  /*//.'(  « - .  h)  -  AV‘(k) 

t  w 

/=— JL* 

(1.22) 

3[  n.  A I  >  U 

J(it.  A)  ^  0 
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where  P  is  a  numerical  dissipation  term.  Since  ^  0.  monotone  .scheme.s  aie 

only  first-order  accurate  approximations  to  the  initial  value  problem  of  Eq  1.9. 

Suppose  the  scheme  given  by  Eq  1.10  is  a  monotone  .scheme  and  thus  provides 
a  second-order  accurate  numerical  approximation  to  the  modified  equation,  Eq  1.21. 
rewritten  as 

+  =  0  (1.23) 

where  g  =  A)ux-  .Applying  this  scheme  to  the  following  equation 


+  (/  +  (l/A)</)j.  =  0 


(1.2d) 


yields  a  second-order  accurate  approximation  to  its  modified  equation.  Since  </  —  0[A.r] 
the  modified  equation  satisfies  [22] 


Ut  -b  A-  =  0  [(A;f)- 


(1.25) 


Thus,  application  of  a  first-order  scheme  to  a  scalar  conservation  law  with  an 
appropriately  modified  flux  function  yields  a  second-order  accurate  approximation 
to  the  original  equation  Ut  +  A  =  0.  Note  that  in  order  to  apply  the  .scheme  to  the 
modified  flux  function,  g  must  be  a  differentiable  function  of  «.  Haiten  achieves  this 
by  smoothing  the  point  values  of  g  (22).  This  smoothing  enlarges  the  suppoit  of 
the  scheme  such  that  his  first-order  scheme,  using  a.  three-point  stencil,  becomes  a 
second-order  scheme  using  a  five-point  stencil.  The  leadei  is  leferied  to  refeience  [22] 
for  the  details  of  how  the  three-point,  first-order  scheme  is  constructed  .so  as  to  ensuic 
its  TVD  property. 

Let  us  now  turn  our  attention  to  the  specific  scalar  scheme  developed  l)y 
Harten.  Consider  a  three-point,  finitc-flifl'etcnce  scheme  in  conservation  form  with 
the  following  numerical  flux  function 

/('^b:  'b+i)  -  ~  ( 1/A)Q  (AA+1/2)  ■2vj+i/2i'|  ( l-‘26) 


9 


where  -  vj  and 


=  . 

Q  is  a  lunction  known  as  the  coefficient  of  luunerical  viscosity.  Numerical  v'isco.sity 
is  the  mechanism  that  allows  a  discontinuity  to  be  captured  as  part  of  the  numerical 
solution  [20).  This  is  in  contrast  to  shock  fitting,  where  the  discontinuity  is  considered 
as  an  internal  boundary. 

Lemma.  3.1  of  reference  [22]  states  that  the  above  scheme  is  TVD  under  the 
Gourant-Friedrichs-Lewy  (CFL)  condition 

AmaN-|rt"^.|/.^|  < /i  (1.28) 

given 

|.T|<Q(.r)<l  (1.29) 

for  0  <  |;!;|  </<<!. 

The  first-order  accurate  three-point  scheme  given  by  Eq  1.26  is  converted  to 
a  second-order  accurate  scheme  b\  appl\iiig  the  three-point  scheme  to  modified  flu.x 
values  J'j‘  (22): 

=  ./■  (vj)  +  ( 1  /A)(/.,  gj  =  fj  (vj-i ,  vj,  Vj^i)  ^  ^ 

~  'G+l/2  +  "6+1/2  ^+1/2  =  (.9^  +  1  ~  .9i) /^J  +  I/2V 

where  i?  =  Xa.  The  modified  numerical  flux  i.s 

■/,;+l/2  ~  2  l-lj'  JjU  ~  ('(i+1/2) 

=  5  [/(^b)  T  ./ (^b+i ))  (L-^O 

-b  (1/(2A))  [</j  +  (jj^x  -  Q  (f^+1/2  +  7J+1/2)  ^j+i/2t''] 
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Lemma  3.2  of  refei'ence  [22]  provides  that  Eq  L.31  represents  tlie  numerical  flu.\ 
of  a  second-order  scheme  so  long  as  Q(.r)  is  Lipschitz  continuous  and  (jj  .satisfic.s 


0}  +  (jj^\  =  ~  ('6+1/2)  j  ^j+i/2^'’  +  O  [2] 


7i+i/2i^i+i/2t'’  =  9}+\  -  <Jj  —  0  [!l‘) 


(1.32) 


Haiien  [22]  constructs  (j  in  the  following  manner  so  as  to  satisfy  Eq  1.32: 


=  .Sj+1/2  ma.x  [0,  min  (|//i+i/2| , ^j-1/2  ■  -6+1/2) ] 
=  Sj-+i/2  min  (j^j+1/2  ,|6_i/2|) 


(^+1/2  ■  6-1/2  >  0) 
(6+1/2  •  6-1/2  ^  0) 


(1-33) 


where 


6+1/2  -  5  |q  ('6+1/2)  ”  ('6+1/2)  j  Aj+i/2y 
■6+1/2  =  6/"  (6+1/2) 


(1.34) 


Finally,  Lemma  3.4  of  reference  [22]  provides  that  a  conservative  finite  differ¬ 
ence  scheme,  with  the  numerical  flu.x  given  by  Eq  1.26,  is  TVD  under  the  lestriction 
of  Eq  1.2s  so  long  as  Q(  t')  satisfies  Eq  1.29.  Thus  a  second-order  accurate  (e.xcept 
near  points  of  e.xtrema  where  Sj+1/2  is  discontinous),  five-point  .scheme  has  been  cuii- 
structed  for  the  .solution  of  Eq  1.9.  The  .scheme  provides  high  resolution  captuiing 
of  discontinuities  and  converges  to  a  physically  relevant  .solution. 

1.3.3  E.riension  to  Syslcnt.^  of  Com-u-rat Ion  Lnwi.  VVe  now  concern  ouiseKe.s 
with  e.xtcnding  the  scalar  scheme  developed  in  .Section  1.3.2  to  systems  of  conserva¬ 
tion  laws.  Currently,  TVD  schemes  are  only  defined  for  scalar  hyperbolic  (onseiva- 
tion  laws  or  constant  coefficient  hyperbolic  .systems.  This  is  due  to  the  fact  that  the 
spatial  total  variation  of  the  solution  to  a  .system  of  nonlinear  conservation  laws  is 
not  necessarily  a  monotonically  decreasing  function  of  time  [46].  Wave  interactions 
may  cause  the  total  variation  to  increa.se.  flarten  extends  the  technique  using  a  gen¬ 
eralized  version  of  Roe's  approximate  Riemann  .solvei  [22].  The  idea  is  to  apply  the 
scheme  in  a  scalar  fashion  to  each  of  the  svstems  linearized  characteristic  variable.^. 


Harten  uses  this  fad  to  extend  his  scalar  scheme  to  general  nonlineai  systems  of 
hyperbolic  conservation  laws. 
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Let  be  the  component  of  =  *j+i  “  '■'j  '’i  the  {R^']  coordinate 

svsteni  such  that 


^i+l/2V  — 

k=l 

Idle  scheme  given  by  Eqs  1.30-1.34  is  extended  to  general  systems  as 

-A 

/i+1/2  =  5(/(^^i)  +  ./’(ui+i)] 

+  ^  T,h=\  -^i+l/2  [(jj  +  t/j+1  ”  Q^'  ('^j+l/2  +  d'j+l/i)  ‘^■;+l/2 
where  i'j^i/2  =  ^<'j+ii2 


with 


Sj 

=  -"j+1/2  (|.^+i/2| ' 3j-i 

1/2  ■  -^'5+1/2)  j 

II 

5  [Q^‘  ('6‘+1/2)  ~  ('6>I/2)  ]  <^J+l/2 

*i+i/2  - 

6+1/2 

(l/i+i  “  3j)  /^i+1/2 

(<^j4-l/2  ^  0) 

= 

0 

(•^1+1/2  =  0) 

(1.43) 


(1.44) 


(1.46) 


(1.47) 


1.3.4  EnU'oinj  Enforcemeid.  .-Vs  a  final  comment  on  the  initial  development 
of  Marten’s  second-order  TVD  scheme,  we  turn  now  to  the  question  of  plnsicallv 
relevant  solutions  for  systems  of  equation.^.  .\.s  mcniioned  in  (he  pievious  .sertiuii. 
the  total  variation  may  not  be  a  monotuiru  ilet leasing  function  of  time  due  to  wave 
interactions.  In  addition.  Oleinik's  entiups  ineijiialiu  ensures  physic all\  iele\aiil.ui 
admissable,  solutions  only  in  the  limit  as  A.c  --  0.  In  realiu  we  are  concerned  uilh 
obtaining  admissable  solutions  on  a  relaii\(4y  coar.se  mesh. 


In  order  lo  arrive  at  a  proper  criterion.  Marten  examines  the  Rieniann  initial 
value  problem  (22j  for  Eq  l.l: 


tt(.i-.0)  =  0(x)  =  ft/,  .r  <  0 

=  u/i  X  >  0 


( 1 .4S) 


with  ft/,  and  ur  satisfying  the  Rankinc-Hugoniot  rclation.s  with  Wd\e  speed  -s.  If 
ft(.i’,  i)  —  o(x  —  $t)  is  to  satisfy  Oleinik’s  inequality  the  numerical  .scheme  mu.st  yield 
a  steady  progre.ssing  profile  with  a  narrow  transition  from  ul  to  hr  [20.  22).  Marten 
refers  to  this  property  as  resoliiUon. 

If  the  solution  u(x,L)  =  (f)(x  —  .s/)  is  inadmissable.  then  the  .solution  is  a  fan  of 
waves  [22].  This  fan  of  waves  is  a  function  of  x/t  and  consi.sts  of  a  rarefaction,  or 
expansion,  wave  in  the  same  field  as  the  initial  discontinuity.  The  plnsical  solution 
requires  the  initial  discontinuity  break  up  instantaeously,  since  u(x,l)  =  0{x/t).Thc 
term  entropy  tnforceincnl  refers  to  the  requiiement  that  the  numerical  scheme  break 
up  the  initial  discontinuity  at  a  fa.il  rate,  thus  imitating  the  physical  behavior  [22]. 

The  systems  of  conservation  laws  under  consideration  contain  two  types  of 
characteristic  fields,  termed  nonlinear  and  linearly  degenerate  by  Marten  [22j.  The 
nonlinear  fields  are  defined  such  that  a^R’^  ^  0,  while  the  lineaily  degenerate  ficid.s 
are  defined  by  =  0.  The  waves  of  a.  nonlinear  field  aie  shock  wa\csoi  expansion 
waves  while  the  waves  of  a  linearly  degenerate  field  are  solely  contact,  or  entropy. 
di.scontinuiiies. 

To  addre.ss  the  que.stion  of  entiopy  enforcement,  (.orisidei  the  .scheme  given  by 
Eq  1.26,  which  has  the  effective  numerical  v'iscosity  coerficicnl 

=  (1-19) 

The  least  dissipative  form  of  0  is  arrived  at  by  choosing  it  to  be  con.si.steni  with 
Eq  1 .29  such  that 

Q{x)  =  |:i] 

With  0  given  by  Eq  1.50,  the  scheme  of  Eq  1.26  can  be  rewriiteu  as  [22] 


jA*-'  =  (,'f  -  (//j+1/2) 


(1.50) 


where 


t/~  =  min(//.  0)  =  ^  (//  -  |//|)  =  max'(7/,  0)  =  3  (/^  -r  |//|)  ( 1 

TLarten  poinlb  out  that  the  scheme  of  Eqs  1.51  and  1.52  is  a  generalisation  of  the 
Couranl,  Isaac.son,  and  Roes  scheme,  which  has  been  thoroughly  analyzed  in  the 
literature.  The  interested  reader  is  referred  to  reference  (22!  for  further  details. 

If  the  scheme  given  by  Eqs  1.51  and  1.52  is  applied  to  the  Riemann  problem 
with  the  Rankine-Hugoniot  relation  satisfied  by  letting  the  speed  of  propagation  be 
zero,  Eq  1.5J  holds  the  initial  discontinuity  steady  regardless  of  cntrop\  considei- 
ations.  In  other  words,  the  initial  di.scontinuity  is  not  broken  up  and  theie  is  no 
entropy  enforcement  in  ihi.s'case. 

The  problem  is  that  the  numerical  viscosity  vanishes  for  //  =  0.  Marten  elimi¬ 
nates  this  problem  by  modifying  Q{x)  =  |.i:|  near  x  =  0  to  be  positive.  The  modifi¬ 
cation  is  as  follows  [22]  for  0  <  c  <  j 

Q(a-)  =  (a:V0k))-}-c  N<2c 

(1.0.3) 

lri:|  |.r|>2e 

with  the  entropy  correction  parameter,  e,  typically  of  order  O.l. 

Marten  summarizes  the  results  of  numerical  experiments  carried  out  with  the 
scheme  of  Eqs  1 . 11  and  1 . 15  apjdicd  to  the  Euler  equations  for  the  Riemann  piobK  in. 
The.se  experiments  itscd  c  =  0.05.  0.1,  and  0.25  for  all  fields,  and  also  c  =■  0  for  ihe 
linearly  degenerate  field.  Basically,  highly  resolved  shocks  were  obtained  foi  all  valiu's 
of  c  under  consideration.  The  contact  surface  was  better  resolved  than  with  the  first- 
order  accurate  scheme  of  Eqs  1.51  and  1.52.  but  still  remained  rathei  smeared. 

To  prevent  excessive  smeaiing  in  the  linearly  degenerate  field  containing  (he 
contact  surface.  Marten  replaces  E(|  l.lfi  in  the  linearly  degeneiate  field  with  [22] 

(1.5-1) 

where  (jj  is  the  right  hand  side  of  t/j  given  by  Eq  1.-16  and.  is 

Tjj  =  S  max  [o,  min  |oj+i/a|)) 


(l.•55) 


II.  Inviscid  Analysis 


2.1  Equations 


The  Euler  equations  are  statements  of  tlie  conservation  laws  for  mass,  momen¬ 
tum,  and  energy  assuming  an  inviscid.  nonconducting  gas.  W'lien  the  Euler  equations 
are  arranged  such  that  p,  pa.  pv,  and  t  are  the  dependent  variables  tiie  conservative 
or  divergence  form  is  obtained.  Lax  showed  that  the  conservative  form  of  the  Eu¬ 
ler  equations  satisfies  the  weak  solution  of  the  Rankine-Hugoniot  relations  and  thus 
correctly  predicts  the  jump  conditions  across  the  shock  discontinuity  [l.  35].  In  fact, 
use  of  the  conservative  form  is  nccessai  v  foi  tlie  discontinuity  to  represent  a  physical 
wave  when  shock  capluring  schemes  are  applit'd  [ij.  The  conserv'ative  form  is  often 
referred  to  as  the  divergence  form  because  the  equations  identify  the  divergence  of 
physical  quantities.  The  governing  equations  may  be  wriu  ..  in  the  following  v'ector 
form: 

dU  ,  dFjU)  ,  dGjU) 

dt  '  dx  ‘  Op 


(2.1) 


where  U  contains  the  dependent  variables  which  arc  the  density,  p\  .i-momentum. 
pu\  //-momentum,  pv.  and  total  energy  per  unit  volume,  e.  F  contains  the  flux  terms 
differentiated  with  respect  to  .x.  and  C  coniains  the  flux  terms  differentiated  with 
respect  to  p.  'J’he  elements  of  V .  F,  and  O  are: 


p 

77? 

11 

in 

/'■’  = 

nr/p  -{-  p 

G  = 

nu 

II 

nn' 

ir/p  -b  p 

e 

(e  -f-  p)iii/p 

(e -b  z?)??//; 

where  tn  =  pii  and  v  —  pv.  I'he  pre.s.siirc.  p.  i.s  given  as 


l>=  (7-  1) 


(nA  -f  »-) 


2p 


(2.3) 


for  a  thermally  and  calorically  perfect  gas. 

.A  general  spatial  transformation  of  the  form  ^  =  ^{.v.p)  and  //  =  pix.p)  i.s  u.sed 
to  transform  Eq  2.1  from  the  phy.siral  domain  {.r.p)  to  the  computational  domain 


{C.,1]).  Tile  sUong  conservation  law  tbrin  of  the  Eulei  e(|iialions  is  now  given  by  [lo 


dU  ,  dF{V)  dG[U) 
dt  '  di] 

(•2.-1) 

0  =  UIJ 

(■2.5) 

F  =  +  iyG)  IJ 

(•2.6) 

C  =  iVrF  +  n,jG)  IJ 

(•2.7) 

J  —  ~  ^y']x 

(2.8) 

where  J  is  the  Jacobian  of  the  transformation. 

Since  the  TVD  method  used  herein  utilizes  the  local-characteristic  approach, 
which  is  a  generalization  of  Roe's  appro.x'imaie  Riemann  .sohcr[36|.  the  .Jacobian.s  .1 
and  B  of  F  and  G  are  required  and  can  be  written  as 

A  =  UrAi-GjB) 

(2.9) 

B  =  (ijrA  -h  >l,jB) 

(2.10) 

wliere 


0  10  0 

„  -^(7  -  1)  {«^  -r  f  ')  -  (3  -  7)1/  ( 1  -  7)r  -  ! 

A  =  /'(,  = 

—  i/f  r  II 

.  [5(7-  !)(»'  +  '•')-//]»  ll~h-\W  ‘.;a 

0  0  10 

—  uv  r  u  0 

B  =  6V  = 

^(7-l)(«*’  +  i’*)-t-'  (l-7)«  (3-7)1’  ';-l 

.  [(h-  +  //-{7-l)r'  7*- 
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with  the  total  enthalpy.  H.  given  by 


The  eigenvalues  and  eigenvectors  of  13  are  obtained  by  replacing  ^  in  Eqs  2.12 
til  rough  2.16  with  //: 


nrU  -r  //;,<•  -  kiC 
-r  Vy^' 

‘{rtl  -T  //„»’  -r  l.-riC 
-r  %{• 


(2.17) 


(2.16) 


(2.20) 

(2.21) 


3.2  .\himeric(il  Procedure 

2.2.1  l-D  Hoc.  I.nx-Wendroff.  mid  TVD  .Mfiorilhms.  The  1-D  scheine.s  under 
consideration  are  the  first  order  accurate  Godunov  -tvpe  [22]  .scheme  of  Roe.  refcni’d 
t.oa.s  the  ROE.scheine:  the  .second-order  acciiiale  h.(.\  W'eiiclrofr-tvpe  .stheme.  referrecl 


Table  2.1.  Dissipation  Terms  for  ROD.  I.VV  and  I'b’iTC  Sclienies 


.Scheme 

Description 

Di.ssip€ation  (.?) 

ROE 

1st  Order  TVD 

‘ 

2nd  Order  non-TV’D 

-  (">-5-1/2)  ^i-fi/2 

ui/nc 

2nd  Order  TVD 

•^il/2  =  0^1/2  *-^1/2)  «5-M/2  -  (^i  -b  UU) 

with  Ecjs  1.54  through  1.5S  applied 
to  the  linearly  degenerate  field 

to  as  the  LW  scheme:  and  the  secoi’.d-ouler  accmaie  TA’D  scheme  of  liarlen,  referred 
to  as  ULTlC.  All  three  schemes  can  be  written  in  the  form 


r;+'  =  v;  -  A 


with  the  appropriate  dis.sipalion  term.  j'.  from  Table  2.1. 


■3.:S.S  3-D  Unrlen-Ycc  Finilc-VoUimc  AlyorUhni.  An  upwind  TV'D  scheme 
in  finite- volume  form  (d.5)  is  used  in  the  pie.'seni  ^tii«ly.  The  grid  spacing  is  denoted 
by  and  A//,  such  that  C  =  and  ;;  =  /.  A;/.  I 't  iiixing  the  .Strang-type  fractional 
step  method  allows  the  scheme  to  be  implemented  in  a  local  characteiislic  approach 
and  ensures  second-order  accuracy: 


where 


#-t/-  .-it  i-it 


rhr'n 

M- 


[•» 


A/ 

A^ 


(2.23) 

{2.21) 


wilii  li  —  A/.  Application  of  the  entire  .sequence  of  operators,  (one  iteiaiion)  advances 
the  solution  two  time  levels.  The  functions  and  are  the  numerical  fluxes 

in  the  ^  and  //  directions  evaluated  at  cell  interfaces.  For  instance.  in  Vee's 

finite- volume  formulation  is  expressed  as 


(/•j.fc  -r  Fj  r\.f:)  -r 


■  M 

‘  M 


(2.26) 


where  the  subscript  j  -i-  ■;  is  a  simplified  noialion  for  j  -r  I’he  numerical  flux 
function  in  the  ;/  direction  is  defined  similarly: 


The  eigenvalues  and  eigenvectors  are  evaluated  at  cell  interfaces  using  sym¬ 
metric  average.s  of  i-jj..  aiul  f  Vt  ^  jj-i-i-  re.spectively.  Roe'.s  averaging 

technique  for  a  perfect  ga.s  i.s  u.se<l  herein  and  Sakt's  the  form  {3Cj 


n  .  I  ..  = 


where 


O-r 

1 

Ov, 

■rl.« 

-i-  I- 

0-r 

1 

Dll 

fT*  l.< 

-F  //,.t 

l>-\ 

-  1 

I 

/  > 

•> 

^  =  yJpj-i-wJPiX 


(2.2S) 

(2.29} 
(2.30i 
1 2.3 1 1 

(2.32» 


Roe’s  averaging  techniqiu*  is  used  bccau.se  it  ha.s  the  computational  advantage  of 
perfectly  rc.solving  .stationary  {dfij,  but  not  necessarih  moving,  disrontiniiities. 


The  cluaiitities  ‘'nd 


formulation 


are  defined  as  follows  for  the  finite-volume 


The  constants  <'iid  necessary  in  determiinng  Eq  2.15,  are 

defined  as 


The  vector  function  is  composed  of  elements  denoted  as  '' 

second-order  upwind  TVD  scheme.  The  elements  are  given  by 


^,4.' 

•  .T+- 


=  <7 


{''m)  W+i  +  ■‘^.0  ■  ^ 


where,  with  A  =  A/./A^ 

0,.  1  is  the  difference  of  the  characteristic  variables  in  the  ^  direction. 

Jf  J 


’  (««  -  bh}/2  ' 

(art  +  bh)/'2 

+ 

•r 

c 

1 

cc 

(2.39) 


(2.30) 
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■where 


■>  I  2 

A^.^ie  + - - -A.^p  -  n.^Aj^im  - 


^»-i-  i.  '  ■ 


cc  =  kiAj^in  +  {koUj^  -  k^v^^)  A^^p  -  k-iA^^m  (2.-13) 

Aj^iz  =  Sj+i,k  -  Zj,k  ('2.44) 

The  difference  of  the  local  characteri.stic  variable.s  in  the  ?/  direction  i.s  obtained  in 
similar  fashion: 

t'n-J  =  Ki  -  '-'m) 

[(<W-«e)/2' 

O'?.,),  A,.,ip-aa 

^+5  _  '>+2'  (2.46) 
{(Id  +  cti)l2 

.<h\  \.  -'f  ■ 

where 


S  -"\+kk  -  %4.i^A-+i'»  -  q-+iAjr.,j.i7)  (2.47) 


2  ^^2 


ee  - - [/.‘i  Ai.,j.im  -  (^ki  Uf.^i  +  k'iVk^i)  +  ^’2A;_.,^in|  (2.48) 

H-+i 

If  =  k\  Af^^in  +  {kMik^i  -  /-i  t’A-+i)  -  k2Ai,,.im  (2.49) 


with 


(^■2)fc+L  = 


The  function  7^  ,  i.  is  given  by 

J'T  T 


where 


■'+2/ 


<r{-v)  =  ^  [Q(x)  -  .r-] 


(2.51) 


(2.52) 


(2.53) 


Q{x)  = 


|x|  (|x|  >  2e) 

{xy{Ac))  +  c  (l:r|<2e) 


The  entropy  correction  parameter,  t,  is  generally  fixed  during  computations,  but  can 
vary  between  0  and  0.5. 

The  function  </'  in  Eq  2.37,  initially  referred  to  in  Section  1.3.2,  i.s  termed  ihe 
‘limiter’  function  and  can  be  exprc.sscd  in  a  variety  of  way.$  [-15].  The  pre.sent  .study 
bases  the  choice  of  the  limiter  on  ihc  t\  pe  of  characteristic  field  under  considc:atic>n. 
For  the  nonlinear  fields,  a//?/  ^  0,  Eq  d.S'ld  of  Yee  [45]  is  used: 


For  the  linearly  degenerate  fields.  a'jV  =  0.  Eq  4.34g  of  Yee  [45]  is  applied: 


f/j  =  S  ■  max  [0,  min  (2  |o'^l|  .  .5’  •  .  min  .  2.S'  ■  o;_,)]  (2.57) 


where 


S'  =  .sgn 


(2.53) 


The  nonlinear  fields  correspond  to  /  =  J  and  /  =  3  while  the  linearly  degeneraic 
fields  correspond  to  /  =  2  and  /  =  4.  It  should  again  be  noted  that  (he  waves  of  a 
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nonlinear  field  are  either  shocks  or  rarefiactioii  waves  while  the  waves  of  a  linearh 
degenerate  field  are  uniquely  contact  discontinuities  [22j.  Since  this  is  a  five-point 
scheme,  the  values  of  are  needed  at  cell  centers  just  outside  the  computational 
domain.  Zeroth-order  extrapolation  is  u.sed  to  obtain  the  necessary  \alucs,  following 
the  example  of  Harten  (22]. 


2.^.2  2-D  Harten- Yec  Chain-Ruk  AUjorithm.  In  addition  to  the  finite- 
volume  formulation  of  Yee,  a  finite-difference  form  based  on  the  chain-rule  con¬ 
servation  form  of  the  governing  equations  w'as  utilized  (-10): 


H  di  ^  ■()(,  '  y;/  ^  di)  ~ 


(2.59) 


Previous  researchers  report  that  the  goxerning  equations  in  this  form  are  more  com¬ 
putationally  efficient  than  the  strong  con.sei\a.lion  form  used  in  the  finite- volume 
approach  [40].  This  w'as  found  not  to  l)e  the  case  for  the  current  TVD  algorithms, 
with  both  formulations  performing  approximately  the  same  in  terms  of  computa¬ 
tional  efiicienev. 


The  local  characteristic  a|)proach  given  by  Eq  2.23  is  now  applied  to  U  instead 


of  U: 


where 


ri»-r-  —  r^'/-  /’ll  /’ll  /’ll  /’IR- 1  :ii 
■  j-l;  ~  '  J.l- 

(2.60) 

(2.61) 

(2.62) 

The  numerical  fluxes,  /'j+r.c.  and  for  ihe  dmin  rule  conservation  form  are 


=  5  +  {(,1,  (6>  +  (2-63) 


26 


— ^/?„  <[),.,  1 

Al  ''••+i  ^■*'2 


The  quantities  (^x)j+i'  (^'i)_,+  L)  and  {k2)j^i  are  denned  as  follows  for  the  chain  rule 
formulation: 


j  +  l.k 


(2.65) 


III.  Inviscid  Results  and  Conclusions 


Chapter  III  details  application  of  the  Marten- Vee  T\’D  algorithms  to  three 
cliircrcnt  classes  of  problems.  Riemann's  problem  of  gas  dynamics  is  covered  first.  .-\ 
shoclv  wave,  rarefaction  wave,  and  contact  surface  are  present  to  test  the  capability 
of  the  TV’D  algorithm  to  resolve  the  features  of  both  linear  and  lineally  degenerate 
fields.  Both  Marten- Yee  algorithms  degenerate  to  Marten's  L’LTlC  scheme  for  this 
case  since  there  are  no  metric  variations.  In  addition  to  the  ULl'lC  .scheme,  solutions 
from  the  Roe  and  La.\-\Vendroff  schemes  arc  presented  to  provide  the  leader  with  a 
performance  comparison. 

.Steady-state  flow  through  a  supersonic  cascade  of  wedges  is  then  c;xamined 
using  the  Marten- Vee  finite-volume  scheme.  Both  shock  and  e.xpansion  waves  are 
-present  in  this  test  case.  Finally,  flow  through  a  typical  transonic  turbine  rotor 
is  considered.  This  test  case  is  used  to  demonstrate  the  capability  of  both  the 
finite- volume  and  chain-rule  algorithms  to  deal  with  transient  stait-up  phenomena 
in  route  to  a  steady-state  solution.  Boundaiy  and  initial  conditions  utilized  for  these 
solutions  aie  di.scussed  at  length  to  procide  an  appreciation  of  how  thc\  complement 
the  physical  nature  of  the  TVD  schemes. 

3.1  Riemann's  Problem. 

Solution  of  Riemann’s  problem  [9]  piovides  ci.  means  for  evaluating  the  ability 
of  a  scheme  to  re.solve  (he  waves  pre.seiit  in  both  nonlinear  and  lincaily  degeneiate 
fields.  The  Roe  (ROE),  Lax'-VVendroff  (lAV),  and  ULTlC  schemes  arc  applied  to 
Riemann's  problem  in  order  to  compaie  Iheii  peiformance.  Solutions  aie  compared 
until  those  of  Marten,  who  utilized  the  same  schemes  [22],  as  a  check  for  correct 
implementation  of  the  algorithms. 


Rieinann’s  problem  is  now  solved  for: 


j  X-  <  0 

\  Un  r  >  0 


(3.1) 


where 


0.4-15 

0.5 

Ul  = 

0.311 

Un  = 

0 

0.8928 

1.4275 

(3.2) 


The.?e  conditions  e.stablish  a  leftward  ino\ing  larefaction  wave,  lightwaid  moving 
contact  surface,  and  rightward  moving  shock  wave.  Figures  3. 1-3.9  show  the  results 
obtained  when  the  ROE,  LW,  and  ULTlC  schemes  are  applied  to  this  problem.  It 
should  again  be  noted  that  ELTlC  is  the  degenerate  form  of  the  Harten-Vee  scheme 
for  the  l-D  problem  with  no  metric  variations.  The  circles  are  the  computed  values 
while  the  solid  line  delineates  the  e.xact  solution.  The  calculations  are  consistent 
with  tho.se  of  Marten  (22)  in  that  they  were  carried  out  to  100  time  steps  with  a  CFL 
restriction  of  0.9.5  using  MO  cells.  In  addition,  a  value  of  e  =  0  was  used  for  the 
TVD  .schemes  with  Roc  averaging  used  only  in  the  ULTlC  scheme.  For  the  one¬ 
dimensional  ca.se.  Roe  averaging  seems  to  be  of  benefit  only  when  acouiaie  re.solution 
of  the  contact  surface  is  desired,  and  then  changes  the  result  only  slightlv  by  biinging 
the  density  at  the  leading  edge  of  the  contact  surface  to  its  correct  value  one  giid 
point  sooner.  X'alucs  of  c  between  0  and  0.2.5  .seem  to  produce  almost  identical 
results  e.xcept  that  c  =  0  seems  to  enhance  the  resolution  of  the  contact  sin  face  in 
the  ULTlC  solution.  This  is  consistent  with  Marten's  observations.  Ocerall.  the 
results  of  this  investigation  seem  to  be  almost  identical  with  those  of  llaiten. 


Figures  3. 1-3.3  show  that  the  first-order  ROE  scheme  provides  a  fair  resolu¬ 
tion  of  the  shock,  but  does  rather  poorly  in  resolving  both  the  rarefaction  wave 
and  the  contact  discontinuity.  Note,  however,  that  the  ROE  scheme  is  T\’D  and 
(hat  its  monotonicity  pre.serving  propertc  prevents  oscillation  of  the  .solution  at  the 
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Figuie  3.1.  Density  from  ROE  Scheme  .Applied  to  Riemann's  Problem 


discontinuities. 

Figures  3. 4-3. 6  show  that  the  performance  of  the  iion-TVD  LVV  scheme  leaves 
much  lo  be  desired.  Not  only  does  the  non-monotonicity  of  the  scheme  cause  se¬ 
vere  oscillations  at  the  contact  and  shock  discontinuities,  but  oscillations  are  also 
occurring  at  the  trailing  edge  of  the  rarefaction  wave,  llarten.  Hyman,  and  Lax  [23] 
point  out  fliat  Lax-VVendrofT  schemes  can  produce  non-physical  solutions,  even  when 
attemi)ts  are  made  to  construct  a  physically  correct  entropy  function. 

Figures  3. 7-3.9  clearly  display  the  improvements  of  the  .second-order  ULTIC 
scheme  over  the  first-order  ffOE  scheme.  The  resolution  of  the  shock,  rarefaction 
waves  and  contact  surface  is  cpiite  good.  It  should  again  be  noted  that  a  value  of 
e  =  0  and  the  use  of  Roe  averaging  are  important  for  resolving  the  contact  surface 
as  accurately  as  possible. 

P’igures  3.10-3.12  show  the  results  obtained  when  ULTIC  is  applied  to  a,  dif¬ 
ferent  set  of  data  for  Riemann's  problem,  the  solid  line  again  representing  the  exact 


THIS 

PAGE 

IS 

MISSING 

IN 


ORIGINAL 

DOCUMENT 


‘1.00 


Figure  3.6.  Pressure  from  LVV  Scheme  Applied  to  Riemaiin’s  Problem 


Figure  3.7.  Density  from  ULTlC  Scheme  .Applied  to  R.iemann’s  Problem 
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Figure  3.10.  Density  from  ULTlC  .Applied  to  the  Ricmaiin  .Shock  Tube 

solution.  This  data  is  ph\'sically  representative  of  a  shock  tube  with 

1  0.12-5 

fJr,=  0  rjn=  0  (3.3) 

2.-5  0.2.50 

The  calculations  in  thi.s  ca.se  were  carried  out  to  -50  time  steps  under  the  CFb  restiii- 
tion  of  0.95  with  100  cells,  consistent  with  Marten  (22).  The  results  show  e.xccllent 
re.solulion  of  the  contact  discontinuiu  a.^  well  as  the  rarefaction  and  .shock  \\a\e.s. 
The  results  appear  to  be  ident  ical  with  i.ho.se  of  Marten. 

S.S  Bouudai'ij  Conditions  for  the  Inviscid  Studies 

Appropriate  boundary  conditions,  in  conjunction  with  initial  conditions  and 
(low  parameters  such  as  Mach  number,  arc  necessary  to  arrive  at  the  particular 
solution  of  interest.  Boundary  condition.s  for  both  the  super.sonic  and  tian.soiiic 
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cascade  flows,  to  be  discussed  in  forthcoming  sections,  are  now  described  in  detail. 


3.:2.l  Inlet  and  Exit  Boundary  Conditions.  If  the  inlet  velocity  is  supersonic, 
all  characteristics  originate  upstream  of  the  computational  bouiidaiy  so  the  four  nec¬ 
essary  flow  quantities  may  be  specified.  Likewise,  if  the  outflow  velocity  is  supersonic 
all  characteristics  originate  inside  the  computational  domain  and  the  four  nece.ssaty 
e.xit  quantities  must  be  e.xtrapolated  from  the  interior.  Second-ortlcr  accurate  e.x- 
trapolation  is  utilised  in  the  schemes  under  consideration. 

Subsonic  inflow  and/or  outflow  presents  a  more  complicated  situation.  In  ap¬ 
plying  the  boundary  conditions  at  the  inlet  and  e.xit  of  the  domain,  it  is  assumed 
that  these  boundaries  are  sufficiently  distant  from  the  Ccascade  so  that  planar  wave 
disturbances  propagate  collinearly  with  the  stream  function.  The  disturbances  arc 
required  to  leave  the  computational  domain  without  reflection,  e.xcept  for  the  re¬ 
flection  of  pressure  disturbances  at  the  e.xit.  For  subsonic  inlet  vcIocitic.s.  the  inlet 
boundary  conditions  are  arrived  at  by  first  assuming  that  the  inlet  is  part  of  an 
imaginary  duct  e.xtending  infinitely  far  up.stream  of  the  cascade.  .All  waves  ladiating 
from  the  comi)utalional  domain  should  pa.ss  the  inlet,  without  reflection,  and  con¬ 
tinue  travelling  upstream  for  all  time.  Specification  of  a  constant  ihermodv namit 
stale  at  upstream  infii'.ity  requires  the  e.xpansion  distmbarice  travelling  upstrtraiM  to 
behave  as  a  simple  wave.  This  I)ehavior  allows  the  application  of  one-dimensional 
characteristic  theory  at  the  inlet  (16). 

For  subsonic  inflow,  only  one  characteristic  runs  from  i  he  inlerioi  of  the  domain 
towards  the  computational  boundary.  Therefore,  three  qiianiil.ie.s  must  be  sfvecified 
while  one  mav-  be  e.xtrapolfitcd  from  the  <loniain  interior.  Fai  upstream,  the  total 
pressure,  and  total  tetnperaturc,  T,^,  are  specified,  while  otily  the  inlet  flow 
angle,  is  specified  at  the  computational  boundary.  The  speed  of  .sound  at  the 
inlet.  C2,  is  extrapolated  from  the  domain  interior.  'Phe  Riemann  itivariant  along  the 


characteristic  spanning  the  expansion  wave  from  leading  to  trailing  edge  is  given  by 


2  .  2 

H - =  V2  + - -C2 

7—1  7—1 


where  V  is  the  magnitude  of  the  velocity  vector.  As  the  velocity  vanishes  far  up¬ 
stream,  the  inlet  velocity  is  obtained  from 


V2  =  (c«>  -  C2)  (3.5) 

7  -  1 

which,  along  with  the  inlet  flow  angle,  determines  u  and  v.  The  inlet  pre.ssure  is 
determined  from  the  isentropic  relation 


Ih  =  p«> 


02  \  27/(7- 1) 


The  speed  of  sound  and  pressure  fix  the  state  point,  uniquely  determining  the  density 
and  internal  energy. 

For  subsonic  axial  Mach  numbers,  simple-wave  theory  is  also  applied  at  the 
exit.  The  exit  is  treated  as  an  open-end  duct  that  exhausts  into  a  plenum,  requiring 
the  exit  pressure  to  match  the  plenum  pre-ssure.  Thus,  all  pressure  disturbances 
are  reflected  back  into  the  computational  domain  from  the  exit.  Two  characteristics 
extend  from  the  interior  of  the  computational  domain  to  the  exit,  while  one  oiiginates 
outside  the  domain.  Thus  only  one  quantity,  in  this  case  pressure,  can  be  specified  at 
the  exii  .\11  other  quantities  must  be  extrapolated  from  the  interioi  of  the  domain. 
The  quantities  chosen  for  extrapolation  are  entropy,  tangential  velocity,  and  the 
Riemann  invariant,  R,nv  The  densitv  is  obtained  from  the  isentropic  relation 


Pz  =  (pzlsint)^^'' 


(3.7) 


where  s,„(  is  the  entropy  extrapolated  from  the  interior.  The  pre.ssure  and  density  fi.x 
the  state  point,  uniquely  determining  the  .si)eed  of  sound  and  internal  energy.  With 


the  tangential  velocity  extrapolated  from  the  interior,  the  axial  velocity  is  obtained 
by  applying  the  Riemann  invariant  in  the  axial  direction: 


=  /?.nw - ^rC3  (3.8) 

7  -  1 

where 

2 

^inv  ~  ^^int  d  Tf'int 

7  -  1 


(3.9) 


and  «„it  and  c„,t  are  the  axial  velocity  and  speed  of  sound  at  the  point  inside  the 
domain  where  the  Riemann  invariant  is  evaluated. 


3.S.S  Periodicity  and  Blade  Boundary  Conditions.  Only  one  blade  passage  of 
an  infinite  cascade  is  analyzed.  Therefore,  periodicity  conditions  are  applied  at  cell 
centers,  or  ghost  points,  located  outside  the  computational  domain.  These  points 
are  lo.cated  along  the  outer  boundary  and  also  along  the  wake  cut  when  a  C-type 
grid  is  utilized.  For  an  H-type  grid,  ghost  points  are  located  along  the  upper  and 
lower  boundaries  upstream  and  downstream  of  the  blade.  At  the  blade  surface,  the 
only  condition  that  can  be  specified  is  the  requirement  for  surface  tangency.  Since 
the  blade  surface  is  mapped  to  a  constant  y  coordinate,  the  normal  component  of 
velocity  is  given  by 


+  VyO 


(3.10) 


while  the  tangential  component  is 


'/ytt  -  iho 

\/'li  + 


(3.11) 


The  requirement  for  surface  tangency  is  met  by  setting 


and 


'4,,.  =  -14,,  (3.13) 

where  j  is  the  ^  index,  0  represents  a  ghost  point  just  inside  the  body,  and  1  is  the 
index  of  the  first  cell  center  above  the  body.  Cell  centers  and  ghost  points  are  used  to 
place  the  blade  surface  along  the  interface  of  the  grid  cell  and  ghost  cell.  This  mesh 
system  helps  ensure  both  consistent  and  conservative  boundary  conditions  [35].  The 
inverse  relation  between  the  Cartesian  velocities  and  Eqs  3.11  and  3.10  then  gives 

Ujfl 
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'hj  >h 

Ko 

\/vl  +  ’fy 

-Vx  Py 

(3.14) 


The  pressure  at  the  ghost  points  is  obtained  by  applying  the  normal-momentum 
equation  at  the  first  line  of  cell  centers  above  the  body  [33]: 

-P  ('/*««  +  VyH)'  =  +  iyVy)  Pi  +  {vl  +  pI)  Pr, 

=  Pn\Jvl  +  vl  (3.15) 

Central  differences  are  used  for  both  the  (,  and  ;/  derivatives. 

One  additional  property  is  needed  to  fix  the  state  of  the  ghost  points.  In  the 
present  study,  an  adiabatic  wall  condition  i.s  chosen  to  provide  this  final  property : 

=0  (3.16) 

■Although  it  is  inconsistent  with  the  Eulei  eqiiatioii.s  to  specilA  eithei  the  tcmpcialuie 
or  its  gradient  at  the  blade  surface  [35].  an  adiabatic  wall  condition  has  been  used 
by  others  [33]  and  yields  results  that  agree  well  with  theory  and  experiment. 
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Figure  3.13.  Grid  for  Cascade  of  Wedges 


S.3  Cascade  of  Wedges 

The  cascade  of  wedges,  previously  examined  by  Denton  [12]  using  his  opposed 
difference  scheme,  is  used  to  demonstrate  the  ability  of  the  finite- volume  T\'D  scheme 
to  capture  well  defined  oblique  shocks.  This  cascade  is  shown  in  Figure  3.13.  The 
cascade  has  an  inlet  Mach  number  of  2.0  and  is  designed  such  that  the  leading-edge 
shock  is  exactly  cancelled  u])on  reflection  to  the  upstream  corner,  resulting  in  uni¬ 
form  flow  between  the  two  parallel  surfaces.  The  grid  used  consists  of  124  points  in 
the  axial  direction  and  40  points  in  the  tangential  direction.  Grid  points  are  clustcied 
near  the  blade  surface  in  order  in  improve  the  accuracy  of  the  solution  to  the  normal 
momentum  equation,  used  to  determine  the  pressure  at  the  wall  (33).  Results  weie 
obtained  using  Ci  =  0.2  for  the  nonlinear  fields  and  0  for  the  linearly  degenerate  fields. 
The  computed  pressure  contours  are  shown  in  Figure  3.14.  The  shock  that  forms 
at  the  leading  edge  is  well  defined  as  is  the  reflected  shock  from  the  lower  surface. 
Cancellation  is  achieved  at  the  upstream  cornei  resulting  in  the  desiicd  uniform  flow 
between  the  parallel  surfaces.  This  is  in  sharp  contrast  to  the  authors'  expciionce 
with  the  .MacCormack  scheme  which,  due  to  shock  smearing,  places  the  reflected 
shock  upstream  of  the  corner  allowing  a  weak  shock  to  be  reflected  back  across  the 
passage  [15].  Figure  3.14  also  shows  the  well  defined  oblique  shocks  occuring  at  the 
trailing  edge  as  a  result  of  the  periodicity  condition.  This  shock  structure  is  simi¬ 
lar  to  the  tramline  shock  structure  (12)  that  can  occur  in  a  transonic  turbine  rotoi. 
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Figure  3.15  contrasts  the  e.xact  and  computed  solutions  in  terms  of  Mach  number 
versus  the  non-dimensional  chord  length.  Circles  and  squares  denote  the  numerical 
solutions  along  the  lower  and  upper  blade  surfaces,  respectively.  The  solid  line  rep¬ 
resents  the  exact  solution.  Unlike  most  second-order  shock  capturing  methods,  the 
TVD  sheme  does  not  display  the  dispersive  errors,  manifested  through  oscillations 
of  the  solution,  that  typically  occur  near  points  where  shock  and  expansion  waves 
are  generated  or  reflected  [15].  An  exact  solution  for  the  expansion  along  the  lower 
surface  was  not  presented  by  Denton,  but  was  computed  by  the  author.  Denton  [12] 
attributes  the  exact  solution  to  Brown  Boveri  L  Co.  of  Baden,  Switzerland. 

3-4  High-Work  Loxu-Aspect-Ralio  Turbine 

The  finite- volume  and  chain-rule  schemes  were  also  applied  to  a  transonic  rotor 
cascade  designed  by  NASA  [42].  The  experimental  turbine  is  a  0.767  scale  model  of 
the  first  stage  of  a  two-stage,  high-pressure  turbine  designed  for  use  in  a  high-bypass 
ratio  engine.  This  model  was  tested  in  the  NASA  Lewis  Research  Center’s  Warm 
Core  Turbine  Test  Facility  [42]. 

Figure  3.16  shows  the  mean-line  velocity  diagram  obtained  from  the  N.ASA 
experiment.  In  the  figure,  V  is  the  velocity  in  a  stationary  frame  of  reference,  lU 
is  the  velocity  in  a  frame  of  reference  moving  with  the  rotor,  and  cr  is  a  condition 
corresponding  to  a  Mach  number  of  unity.  Subscripts  1,  2,  and  3  correspond  to 
the  stator  exit,  rotor  inlet  and  rotor  exit  respectively.  Using  the  mean-line  blade 
coordinates  from  reference  [42],  and  the  relative  gas  angles  fiom  Figure  3.16.  the  C- 
type  grid  shown  in  Figure  3.17  is  coiisti noted,  d'hegiid  is  made  up  of  177x20  points, 
again  with  points  clustered  near  the  surface  foi  impioxed  accuraev  in  apphing  the 
boundary  conditions.  The  rounded  trailing  edge  is  replaced  by  a  cusp  to  i^reveiit 
the  severe  expansion  around  the  blunt  configuration.  This  was  found  not  to  be  a 
necessary  requirement  on  the  trailing  edge  geometiy,  and  solutions  were  obtained 
for  a  grid  where  the  rounded  trailing  edge  was  left  intact.  Points  are  also  clustered 


at  the  leading  and  trailing  edges  for  improved  resolution.  21  points  are  placed  along 
the  portion  of  the  C-type  grid  representing  the  inlet. 

The  inlet  and  exit  condions  necessary  for  input  to  the  code  are  derived  from 
the  NASA  test  data  [42]  given  in  Table  3.1.  In  the  table,  t  identifies  total  properties 
and  R  denotes  a  frame  of  reference  moving  with  the  rotor.  Conditions  at  station  1, 
the  stator  exit,  are  taken  to  be  the  conditions  also  existing  at  station  2.  the  rotor 
inlet. 

Consistent  with  subsonic  inflow  at  the  computational  inlet,  the  total  pressure 
and  total  temperatifre  in  the  quic-scent  region  infinitely  far  upstream  of  the  cascade 
are  required  as  boundary  conditions.  The  values  used  are  =  29.32  x  lO'CV/n?" 
and  Tt^  =  420. 2A'.  The  static  pressure  at  station  3,  the  rotor  exit,  is  input  as  the 
exit  pressure.  In  particular,  =  12.05  x  lO'A’/m^. 

3.4.1  Numerical  Soluliou.  The  initial  conditions  applied  for  the  present  study 
are  referred  to  as  “cascade  tunnel  start"  conditions  because  of  the  analogy  to  the 
starting  of  a  blow-down  cascade  tunnel.  The  domain  is  initialized  at  zero  velocity, 
the  pressure  and  temperatuie  corresponding  to  that  in  the  quiescent  region  upstieam 
of  the  inlet.  This  is  analogous  to  placing  a  diaphragm  at  the  exit  of  the  computa¬ 
tional  domain.  At  time  to,  the  .solution  is  started  and  a  centered  expansion  wa\e 
propagates  upstream.  It  is  also  po.ssible  to  place  the  diaphragm  anywhere  in  the 
computational  domain,  but  placing  it  at  the  exit  avoids  the  formation  of  a  contact 
surface  that  must  pa.ss  through  the  domain.  While  the  present  T\’D  scheme  has 
demonstrated  the  ability  to  re.solve  siuh  ,1  <  out  act  surface  in  very  fine  detail,  con¬ 
vergence  is  slowed  due  to  the  fact  that  iIk'  contact  surface  progre.sses  through  the 
domain  at  the  convective  velocity. 

When  the  cascade  tunnel  start  is  u.sed  and  the  diaphragm  is  placed  at  the  exit  of 
the  computational  domain,  a  centered  (Wpaiision  \va\e  propagates  upsteam  through 
the  blade  pa.ssage  and  towards  the  inlet.  .As  the  leading  edge  of  the  expansion 


Table  3.1.  N.^S.A  Turbine  Test  Data 


7  =  1.4 

Ratio  of  Specific  Heats 

Tt,  =  422.2  K 

Inlet  Total  Temperature  (Absolute) 

Pi,  =  31.03  X  rn'N/m'’- 

Inlet  Total  Pressure  (.Absolute) 

II 

<P 

Inlet  Total  to  Static  Pressure  Ratio  (Absolute) 

Phn/Pi  =  1-652 

Inlet  Total  to  Exit  Static  Pressure  Ratio  (Relative) 

Ph/Ph  =  ■-■•160 

Inlet  Total  io  Exit  Total  Pressure  Ratio  (Absolute) 

{V/Vcr)o  =  0.88S 

Inlet  Critical  Velocity  Ratio  (.Absolute) 

{V71/,)3  =  0.384 

Exit  Critical  Velocity  Ratio  (.Absolute) 

(H7HCr)2  =  0.381 

Inlet  Critical  Velocity  Ratio  (Relative) 

(H7H/,,)3  =  0.841 


E:<il  Critical  Velocity  Ratio  (Relative) 


wave  reaches  the  leading  edge  of  the  airfoil,  circulation  is  established  around  the 
blade  through  the  shedding  of  a  starting  vortex  from  the  airfoil.  Since  vorticity  is 
related  to  the  entopy  gradient  per  Crocco’s  equation,  entropy  contours  can  be  used 
to  highlight  regions  of  vorticity.  Figure  3.18,  a  plot  of  the  entropy  contours  after 
1000  time  steps,  clearly  shows  the  starting  vortex  that  has  been  shed  by  the  airfoil. 
The  vortex  is  convected  downstream  and  eventualK  exits  the  computational  domain 
without  being  reflected.  Figure  3.18  also  provides  a  graphic  representation  of  the 
periodic  behavior  in  a  cascade  flow. 

Steady-state  solutions  obtained  with  the  finite- volume  TVD  formulation  com¬ 
pares  extremely  well  against  the  experimental  mean -line  data.  Computational  data 
is  given  in  Table  3.2  for  a  C-type  grid  with  357  x  40  points.  Grid  spacing  is  roughly 
half  that  of  the  grid  shown  in  Figure  6.42  and  is  utilized  in  an  effort  to  verify  the 
location  of  the  stagnation  point  at  the  leading  edge.  Since  the  stagnation  point  is 
determined  by  the  circulation  around  the  airfoil,  the  entire  grid  was  refined  rather 
than  just  the  area  in  the  vicinity  of  the  stagnation  region.  Computed  values  are 
compared  against  experimental  values  through  the  percent  difference: 

%Diff  =  -  1^  X  100  (3.17) 

When  the  .solution  process  was  performed  using  the  chain-rule  formulation 
no  noticeable  difference  in  the  computed  quantities  was  observed.  This  level  of 
agreement  provides  an  excellent  argument  for  the  use  of  T\’D  schemes  in  computing 
transonic  cascade  flows. 

3.5  Conclusions  Based  on  fnviscid  InoesUfjations 

TVD  schemes,  because  of  their  foundation  in  the  mathematics  and  physics  of 
hyperbolic  conservation  laws,  are  clearlv  .superior  to  the  widely  used  Lax-V\’endroff 
family  of  schemes  that  solve  the  partial  diffeicntial  equations  without  regard  to 
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.  Entropy  Contours  Showing  Starting  Vortex 


Table  3.2.  Computed  Versus  Experimental  Data 


,  Quantity 

Experiment 

Computed 

%  Difference 

Pi2rIP^ 

1.652 

1.643 

-0.-54 

Phn 

18.64  X  10‘DV/m2 

18.98  X  10“AVm2 

1.82 

(HVH'cr)., 

0.381 

0..383 

0..52 

(H7M'cr)3  ■ 

0.841 

0.853 

1.43 

Exit  Angle,  6z 

-66.-50° 

-66.75° 

0..3S 

selecting  the  physically  meaningful  particular  solution.  The  superior  performance  of 
TVD  schemes  in  resolving  rarefaction  waves,  contact  discontinuities,  and  shock  waxes 
was  demonstrated  using  Riemann’s  problem.  Solutions  for  both  super.sonic  and 
transonic  flows  exhibit  gieatly  improved  resolution  over  second-ordei  La.\-\\’endiuff 
type  schemes  used  previously  [Ld,  16).  Solutions  also  compare  favorabh  with  the 
available  experimental  and  analytical  data. 

One  important  improvement  not  mentioned  previously  is  that  the  down.bt  ream 
periodic  boundaries  for  the  turbine  analysis  do  not  have  to  be  treated  a.s  solid  walls 
at  any  point  during  the  solution  process.  Previous  experience  of  the  auLlioi  and 
others  (15,  16,  .38)  showed  that  numerical  difficulties  are  encountered  with  the  .Mat- 
Cormack  scheme  if  the  cascade  tunnel  start  is  used.  The  downstream  boundaries  had 
to  be  treated  as  solid  w'alls  until  the  solution  evolved  to  a  point  whcic  the  flow  i)e 
came  aligned  with  the  channel.  .No  such  difiiculty  has  been  observed  with  eithei  tin- 
finite- volume  or  finite-difference  TVD  schemes  described  in  Sections  2.2.2  and  2.2.3. 


One  observation  should  be  made  regarding  grid  cell  skewness.  While  no  adveibC 
effects  due  to  varying  aspect  ratio  were  observed  for  either  the  finite- volume  or  finite- 
difference  formulations,  excessive  cell  skewness  in  the  turbine  cascade  grid  led  to  a 
rather  high  degree  of  entropy  production.  This  appears  to  be  unrelated  to  boundaiy 
conditions  since  the  production  was  most  noticeable  in  the  interior  of  the  domain 
above  the  suction  surface,  where  the  grid  tended  to  be  most  skewed.  Reducing  t 
to  0  tended  to  alleviate  the  problem,  suggesting  that  increased  c  values  tend  to 
magnify  the  effect  of  numerical  viscosity  generated  by  cell  skewness.  The  problem 
was  observed  with  both  the  finite- volume  and  finite-difference  formulation.s,  although 
the  finite-volume  formulation  tended  to  enhance  the  production  of  entropy.  When 
the  skewness  was  reduced,  no  noticeable  differences  in  the  solutions  using  either 
formulation  were  observed.  No  significant  variation  in  the  solutions  is  observed  for  t 
values  ranging  from  0  to  0.4  in  the  nonlinear  fields,  so  long  as  skewness  is  kept  to  a 
minimum.  A  value  of  c  =  0  was  consistently  used  for  the  linearly  degenerate  fields. 


CFL  numbers  as  high  as  0.95  were  consistently  used  to  obtain  steady-state 
results.  In  fact,  the  CFL  number  was  dropijed  to  0.5  only  if  a  contact  surface  was 
in  the  vicinity  of  the  rounded  trailing  edge  of  the  blade.  .At  all  other  times  the  C-FL 
number  was  maintained  at  0.95.  This  is  in  contrast  to  CFL  numbers  as  low  as  0.2 
required  during  startup  and  onh  as  higii  as  O.S  to  maintain  stability  when  using  the 
MacCormack  scheme  [15,  16]. 


The  data  processing  rate  is  1.2425  x  10"'’  .seconds  per  grid  point  per  time  level 
for  the  chain-rule  formulation,  and  1.227 1  x  10"^  .seconds  per  grid  point  per  time  level 
for  the  finite-volume  formulation:  the  piote.s.sing  rate  refers  to  the  C'R.\5'  .N-.\IP/216 
computer.  The  solution  is  monitored  until  laU  ulations  consistently  .-.how  less  than  a 
0.02%  change  in  the  total  energy.  The  lime  dependent  solution  is  then  considered  to 
have  asymptoted  to  the  steady-state  .solution.  .A  typical,  time-accurate  calculation 
requires  approximately  4000  oircrator  .sweeps  to  achieve  steady  slate  convergence. 
This  translates  to  approximately  o.t'  minutes  of  CR.AV  X-MP/216  CPF  time  foi  the 
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177  X  20  grid.  When  a  local  time  stepping  procedure  is  used,  appro.ximately  2000 
time  steps  are  required  to  achieve  the  same  level  of  convergence.  Thus  the  CPU  time 
is  reduced  by  a  factor  of  two.  .A  description  of  the  .ATEC  routines  is  presented  in 
.Appendix  B.  .Appendix  B  also  summarizes  the  rc.sults  of  the  CR.AY  FLOWTR.ACE 
option  used  to  obtain  a  relative  performance  evaluation  of  the  routines. 


IV.  Introduction  to  Part  II 


/,.!  TVD  Schemes  Qiid  the  Navier-Stokes  Equations 

Soon  after  the  in  .loduction  of  the  TVD  methodology  by  Harten  [22]  the  scheme 
began  to  be  applied  to  the  Navier-Stokes  equations.  The  earliest  application  known 
to  the  author  was  by  Chakravarthy  et  al.  [7]  in  1985,  followed  by  Miiller  [.32]  in 
1989,  Riedelbauch  and  Brenner  [34]  in  1990,  and  Lin  and  Chieng  [25],  Seider  and 
Hanel  [39],  ?rd  .Josyula,  Gaitonde,  and  Shang  [24]  all  in  1991.  .All  these  investigations 
dealt  e.xclusively  with  the  steady-state  problem.  Numerous  other  researchers  have 
undoubtedly  applied  the  TVD  methodology  to  the  Na.vier-Stoke.s  equations,  but  the 
author  is  mainly  aware  of  the  above  efforts.  Most  of  the  effort  has  been  directed 
toward  the  investigation  of  hypersonic  flows,  but  Lin  and  Chieng  and  Seider  and 
Hanel  have  investigated  the  transonic  regime  through  solutions  of  the  thin-layer 
Navier-Stokes  equations. 

The  present  effort  is  an  attempt  to  extend  the  application  of  the  TVD  method¬ 
ology  in  two  directions.  The  first  direction  is  the  calculation  of  unsteady  flows,  where 
the  time  accuracy  of  the  scheme  is  important.  The  second  direction  concerns  flows 
wdiere  comi)lex  wave  pheomena  are  present,  but  are  relatively  weak  compared  to 
those  of  previous  investigations.  A  primary  assumption  of  previous  investigations 
is  that  the  flows  are  dominated  by  inviscid  effects:  moderate  or  strong  shock  wares 
arc  pre.sent  in  the  flowfield.  This  assumption  allowed  investigators  to  conclude  that 
solutions  far  away  from  the  boundary-layer  are  accurate.  c\eii  though  the  effect  of 
the  TVD  dissipation  terms  on  the  true  viscositj  in  the  boundary  layer  rernairred  un¬ 
known  [45].  Seider  and  Hanel  w'ere  the  first  to  investigate  the  effect  of  this  dissipatiorr 
on  the  boundary  layer  and  the  present  work  attemps  to  extend  this  kitowdedgc. 
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Jf.S  Overview  of  Part  II 


The  present  effort  is  primarily  concerned  with  the  development  of  an  algorithm 
capable  of  analyzing  laminar  flows  with  boundary-layer  separation  and  heat  transfer 
induced  by  both  steady  and  unsteady  shock  waves.  Several  algorithms  have  been 
developed  that  are  reasonably  accurate  in  predicting  pressure  distributions  for  the 
laminar  shock-boundary-layer  interaction  problem.  To  the  author's  knowledge  no  al¬ 
gorithms,  other  than  those  developed  herein,  currently  exist  that  accurately  predict 
skin  friction  coefficients  in  the  interaction  region  or  the  correct  .sepaiation  and  reat¬ 
tachment  locations.  Similarly,  no  algorithm  is  available  that  accuiately  computes 
local  heat  flux  for  even  the  simplest  geometries  when  shock  wa\es  impinge  upon  the 
boundary  layer. 

To  accurately  compute  the  complex  flow  structure  of  shock-induced  boundary- 
layer  separation,  or  compute  accurate  heat  flux  levels,  the  algorithm  must  provide 
for  high  resolution  of  the  complex  wave  systems  and  maintain  the  proper  physical 
behavior  of  the  problem  under  consideration.  TVD  schemes,  which  lend  themselves 
to  limited,  but  extremely  rigorous,  analysis  provide  the  best  foundation  to  build 
upon.  Although  developed  for  the  solution  of  scalar  hyperbolic  conservation  laws. 
TVD  schemes  perform  well  on  systems  of  hyperbolic  equations,  such  as  the  Riemann 
problem  analyzed  in  Part  I.  The  T\'D  methodology  is  adapted  herein  to  provide 
accurate  .solutions  to  the  parabolic  Navier-Stokes  equations. 


Part  11  begins  with  the  casting  of  the  Navier- .Stokes  equations  in  conservative 
form,  after  which  the  svstem  is  linearized.  'l'\.o  versions  of  TVD  algorithms,  the  1st- 
Order  AFIT  TVD  Navier-Stokes  Code  (A'l'NSC'l)  and  the  2nd-0rder  A  FIT  TVD 
Navier-Stokes  Code  (ATNSC2).  are  then  developed.  Both  algorithms  aie  extensions 
of  the  Harten-Yee  inviscid  algorithm  outlined  in  Section  2.2.3.  .A.TNSC1  is  formally 
first-order  accurate  in  time,  second-oider  accurate  in  space.  .ATNSC2  is  formally 
second-order  accurate  in  time  and  space.  .-VTNSCi  and  .ATNSC2  arc  first  applied, 
along  with  a  fiax-VVendroff  algorithm,  to  the  viscous  Burgers'  equation  as  a  test  ca.se. 


This  test  case  illustrates  the  superior  performance  of  the  ATNSC  schemes,  as  well  as 
the  necessity  of  utilizing  the  fully  second-order  ATNSC2  algorithm  for  low  Reynolds 
number  flows.  The  ATNSC  algorithms  are  then  applied  to  the  solution  of  the  shock¬ 
boundary-layer  interaction  problem.  Computed  solutions  are  compared  with  the 
experimental  data  of  liakkinen  et  ah,  and  with  solutions  obtained  from  Visbal’s 
Beam- Warming  algorithm  [47],  in  order  to  illustrate  the  superior  performance  of  the 
TVD  based  algorithms. 

The  ATNSC  algorithms  are  next  applied  to  the  problem  of  unsteady  shock- 
induced  heat  transfer.  Solutions  are  compared  with  those  obtained  from  VisbaTs 
Beam- Warming  algorithm,  the  iheoiy  of  Mirels  [30],  and  the  experimental  data  of 
Smith  [41].  ATNSC  solutions  are  shown  to  behave  in  a  physically  correct  manner, 
providing  extremely  accurate  solutions.  The  Beam- Warming  algorithm  allows  the 
formation  of  nonphysical  waves,  including  expansion  shocks,  for  this  test  case. 

Finally,  conclusions  arrived  at  from  the  current  investigation  are  summarized. 
Suggestions  are  also  given  for  further  research  involving  use  of  the  ATNSC  algo¬ 
rithms. 
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V.  Viscous  Analysis 


5. 1  Navier-Stokes  Equations 

The  conservative  form  of  the  Navier-Stokes  equations  is  written  as 

du  dFjU)  dGjU)  _  dF.iU,U,,Uy)  dG,iU,U,,Uy) 

dt  dx  dy  dx  dy 

where  U,  F,  and  G  are  the  same  as  for  the  Euler  equations,  Eq  2.1.  and  are 
the  viscous  flux  terms,  given  as 


-|-  VT^y  —  q^ 


UTxy  +  'OTyy  -  qy 


Tary,  and  Tyy  are  the  viscous  stresses: 

Txx  =  (2/i  -h  \)Ux  +  XVy 

Txy  —  {J-  (tiy  -}*  Ui)  (^‘3) 

Tyy  -  (2/i  4-  A)t;y  "f  AUj 

where  y  and  A  are  the  first  and  second  coefficients  of  viscosity  respectively.  The  first 
coefficient  of  viscosity  is  determined  using  Sutherland’s  formula  [1]; 


11  =  Cy 


T  +  C2 


where  Cy  =  1.458  X 10  ®  kg/  (m  -  s  •  and  C2  =  110.4  K.  The  second  coefficient 
of  viscosity  is  given  by 

B  =  2  +  -  (5.5) 
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where  B  =  4/3  yields  Stoke’s  hypothesis,  A  =  —2/3/i.  Solutions  were  also  arrived 
at  using  5  =  2,  based  on  Sherman’s  work  as  reported  by  White  [49].  No  difference 
Wcis  observed  in  the  numerical  solutions  using  5  =  4/3  or  5  =  2. 

The  quantities  qx  and  qy  are  components  of  the  heat  flux  vector,  q  =  —kVT. 
The  coefficient  of  thermal  conductivity,  k,  is  determined  from  the  Prandtl  number, 
Pr: 

Pr=^  (5.6) 

k 

with  Pr  =  0.72  for  air. 

The  equations  may  be  written  in  linearized  form  as 


Ut  +  AUx  +  BUy  =  A\Ux  A-  B\Uy  A-  A^Uxx  +  B2Uyy  A  (As  +  5$)  Uxy  (5’"^) 


where  the  viscous  Jacobian  matrices  are 


Ai  =  dFJdU  A2  =  dFyldUx  A3  =  dFJdUy 
Bi  =  dG^ldU  B2  =  dGyldUy  B3=^dGJdUx 


(5.8) 


with  the  individual  terms  given  in  Appendix  A. 

A  general  spatial  transformation  of  the  form  ^  =  ^(x,  y)  and  y  =  y{x,  y)  is  used 
to  transform  Eq  5.7  from  the  physical  domain  (x,  y)  to  the  computational  domain 


Ut  +  AU^  +  55,,  =  AiU^  A  B\Uy  A  +  525,,,,  +  ^As  +  53^  5^,,  (5.9) 


where 

A  =  ^lA  +  iyB 

Al  =  ifxAl  +  iyB\ 

A-i  =  (iA2  +  CjB2  +  ^x^y  (As  +  Bz) 

As  ~  ^xVyAz  A  ^ijVxBz  A  ^xVxAz  A  ^yVyUz 


(5.10) 
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and 


t/iA  +  rjyB 
VxAl  +  TjyBl 

vIA2  +  VIB2  +  VxVy  (As  +  B3) 
(xVyBz  +  ^yVxAs  +  ^xVxA2  +  ^yVyB2 


(5.11) 


5.S  Numerical  Procedure 


5.2.1  Ist-Order  Time,  2nd-0rder  Space  Algorithm.  A  first-order  time,  second- 
order  space,  upwind  TVD  scheme  is  now  presented  for  the  Navier-Stokes  equations. 
Based  upon  the  excellent  results  achieved  in  the  inviscid  case,  a  chain-rule  formula¬ 
tion  is  utilized.  The  scheme  for  the  Euler  equations  ,  described  in  Section  2.2.3,  is 
second-order  accurate  in  space  and  time.  Taylor  series  expansion  shows  the  scheme 
is  a  representation  of 


Ut  -|-  +  VxFji  +  ^yG^  -1-  rjyG,^  =  ^  —Un  -f  A^U^^  -f  (AjB  -1-  jBA)  U^ri 


and  is  second-order  accurate  for  the  Euler  equations,  since 


u„  =  Ayjii  +  [Ab + bA)  Uf„  +  B% 


(6.12) 


(5.13) 


Viscous  terms  are  added  to  the  Euler  scheme,  Eqs  2.61  and  2.62,  using  second- 
order  accurate,  central-difference  approximations: 


=  Vlu  -  Pj.,)  -  A(tJ 


(5.14) 


fen  - 


(5.13) 
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where  F  and  G  are  given  by  Eqs  2.63  and  2.64.  The  viscous  terms,  and  are 


K«,.  -  V...)  +  (g.,*,..  -  G., (5'16) 


and 


The  scheme  given  by  Eqs  5.14  through  5.17  is  a  representation  of 


Ut  +  ixF^  +  TJxF^  +  +  TjyGy 


E'Xamination  of  Eq  5.9  reveals 


—  ^xFv^  +  VxFvy  +  ^yGy^  +  TjyGy^ 

+f  [-Uu  +  A^U^i  +  [aB  +  BA) 

■  +PUr,r,]+0[AF,Ae,^n^] 

(5.18) 


Uu  =  {P  +  B^,-m-BiB)U,, 

+  ^A\Bi  +  BiAi  •—  AiB  —  BAi  —  i4J5i  —  ByA  +  AB  +  B/^  U^r) 

+  {B\P  +  pBl  —  Bp  —  ppj  UrjtjT,  +  pUrjT^rtr) 

+  {^BzBy  +  Byp  —  pB  —  Bp  +  ByAz-\-  pBy 

—BAz  —  AzB  +  BzAy  +  Ayp  —  pA  —  AB'^ 

+  [BzP  +  PP  +  PP  +  PP'j  U^rjyry 
+  {B^  +  /I3  +  A2P  +  PP  +  PP  +  53^3) 

/aa  Aa  aa  aa  aa  aa 

+  (/IsAi  +  AyAz  ~  .43/I  —  /I/I3  +  /1ijB3  +  BzAy 
— ABz  ~  pA  +  AzBy  +  ByAz  ~  AzB  —  BA-^ 

+  ^v42.43  +  .43i42  +  AzP  +  ^3.42^  Pi^ri 

+  (a^  +  a?  -  AA,  -  A,a) 

+  ^AiA2  +  A2A1  —  AA2  —  A2A)  (5.19) 
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Since  the  term  of  0[Ai]  in  the  truncation  error  of  Eq  5.18  does  not  vanish  upon 
substitution  of  Eq  5.19,  the  scheme  obtained  by  adding  central-difference  represen¬ 
tations  of  the  viscous  terms  to  the  Euler  scheme,  Eqs  5.14  through  5.17,  is  first-order 
accurate  in  time  and  second-order  accurate  in  space.  This  new  scheme  represents 


Ut  +  -1-  r)xFr,  -i-  iyGi  -b  r]yGrt  =  +  rjyGv,,  ,,  , 

(o.iUj 

This  first-order  time,  second-order  space  scheme  is  hereafter  referred  to  as  ATNSCl 
which  is  shorthand  for  Ist-Order  AFIT  TVD  Navier-Stokes  Code.  A  description 
of  the  ATNSCl  routines  is  contained  in  Appendix  C.  Also  included  is  a  relative 
performance  evaluation  of  the  routines,  obtained  using  the  CRAY  FLOWTRACE 
option. 


5.S.2  2nd-0rder  Time,  2nd-0rder  Space  Algorithm.  All  known  TVD  solu¬ 
tions  to  the  Navier-Stokes  equations,  prior  to  the  present  effort,  have  implemented 
the  viscous  terms  in  a  manner  analogous  to  that  of  Section  5.2.1.  A  second-order 
accurate,  upwind  TVD  scheme  is  now  developed  for  the  Navier-Stokes  equations.  As 
previously  mentioned,  the  scheme  for  the  Euler  equations  is  second-order  accurate 
in  space  and  time  and  is  a  representation  of  Eq  5.12.  Utilizing  the  fractional  step 
method,  consider  a  scheme  of  the  following  form  for  solution  of  the  Navier-Stokes 
equations: 


(6.21) 


where 


=  (l  -  AtASi  +  ^AH^]  C/?,  +  A(-P(” 


(5.22) 
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and 


rhTTn 


=  (l  -  AiB«,  +  Ul,  +  A»; 


(5.23) 


where  h  =  At.  The  numerical  viscous  flux  derivative,  is  the  representation  of  the 
viscous  derivatives  terms  plus  terms  necessary  to  cancel  any  first-order  truncation 
error.  Above,  Sf^  represents  a  second-order  accurate  (centered)  difference  approxi¬ 
mation  of  the  derivative  with  respect  to  /.  The  functions  and  are 

the  same  numerical  fluxes  as  for  the  Euler  equations,  Eqs  2.63  and  2.64. 

CfU  provides  a  second-order  accurate  solution  to  the  one-dimensional  equation 


U,  +  ic/f  =  -  yC/,,  +  +  0  [Ai^  Af 

while  C!^U  provides  a  second-order  accurate  solution  to 

Ut  -h  BUr,  =  -1-  -}-  O  [A^^  At?^] 


(5.24) 


(5.25) 


In  the  two  equations  above,  ^  is  the  exact  representation  of  the  viscous  flux  deriva¬ 
tives  plus  terms  necessary  to  cancel  any  first-order  truncation  error. 

Subtracting  selected  viscous  terms  from  the  two  equations  above  gives  quasi- 
one-dimensional  forms  of  the  Navier-Stokes  equations  on  the  left-hand  sides: 


Ut  +  AU^-AiUi-A^U^i-AzU^,  =  ^^-fUu  +  fAW^^ 

— AiC/^  —  A^U^^  —  AzU^n 


(5.26) 


BUr,-  BxUr,-  BzU,,-  BzU^n  =  -  f  t/u  +  f 

—  B\Ur,  —  BzUr,,,  —  BzU^r, 


(5.27) 
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The  objective  is  to  select  ij)  such  that  the  right-hand  sides  of  Eqs  5.26  and  5.27 
equal  zero,  yielding  the  quasi-one-dimensional  Navier-Stokes  equations: 


Ut  -f  AO^  —  AiU^  —  A2U^^  —  AzU^jj  =  0  (5.28) 


Ut  4-  BU,  -  BiU„  -  B^U,,  -  B^U^r,  =  0  (5.29) 

Examining  the  quasi-one-dimensional  Navier-Stokes  equation  in  the  ^-direction  gives 

Uii  =  —  AAi  —  Ai  A  A^ 

-b  ^Aii42  -{-  A2A1  —  AA2  —  i42.i4^  "h 
-f  -{-  A^Ai  —  AA3  — 

+  (A2A3  -b  A^A^  '  (5.30) 

Substitution  of  Uu  into  Eq  5.26,  and  setting  the  left-  hand  side  equal  to  zero,  yields 

+  Y  ((^1  “  ~  ^1^) 

-b  (A1A2  -b  A2A1  —  AA2  —  A2/i) 

+  (A1A3  -b  A3/I1  —  A  A3  —  A3  a) 

+  (A2A3  -b  A3A2)  (5.31) 

Similar  manipulation  of  the  ?7-direction  equation  gives 

^  +  Y  [(jsf  -  BS.  -  b.b) 

+  {B\B2  -b  B2B1  —  BB2  —  B2B^  C7ijt}ij  + 

+  i^BiBs  +  B3B1  —  BBz  —  BzB^  U^,jri 

-b  (44  +  44)  (5.32) 
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Applying  the  fractional  step  operators,  as  in  Eq  5.21,  with  the  numerical  vis¬ 
cous  flux  derivatives  given  by 

+  (Aii42  +  •A2A1  —  i4A2  ~  ^2^4^  -b 
+  (AiAs  +  MM  —  AM  —  AsA) 

+  (A2A3  +  A3A2)  U  (5.33) 

and 

-{-  [BxB2  +  ^2^1  -  BB2  -  B2B)  5l  -{-  Bl6‘l 
+  (.^1.^3  +  B3B1  -  BB3  -  B^B^ 

+  {B2B3  +  B3B2)  5|„„„  +Bl6^^^^]  U  (5.34) 

results  in 


f'h  fh  rh  rh-rin 


{  1  -i-  2Af  {  (A  -  5)  5,  -b  [(b^  -b  Bl  -  BB,  -  B,B)  At  -b  A 
+  [A  +  A  +  (AiA  +  A  Ai  -  AiB  -  BAx  -  ABi  -  A  A  -b  AB  -b  BA)  At)  5, 

-b  (A  A  +  A  A  —  BB2  —  A^)  At5^  -b  B\At8‘l^ 

+  (A  A  +  A  A  —  MB  —  BBz  +  A  A3  +  A3  A 
—BAz  —  A3B  -b  AAi  +  Ai  A  “  Aa  —  A  A)  At^^, 


-b  (A  A  +  A  A  +  Aas  -b  A3  A)  At8'^^^^ 

d*  (A  d"  A3  -b  A2A  +  MM  +  MM  "b  MA^  At8'^^^jj 
-b  (A3A1  +  A1A3  —  A3  A  —  A/I3  -b  /ii  A  "b  AAi 

—AM  —  A  A  -b  A  A  d-  A  A2  —  A2B  —  BA-^  At8^^jj 
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+  {^A^Az  +  •<43^2  +  AzBz  +  BzA-^ 

+  +  [(■^^  4"  ^^1  “■  AAi  —  Ai/lj  Ai  +  -421 

+  ^AiA2  +  4.24,1  —  442  —  424^  Af5|  +  42Ai5|^J’  Ulk 
■hO[At^,A^\Av^)  (5.35) 


The  above  scheme  is  a  second-order  accurate  representation  of 

Ut  -f  AU^  -1-  BUjj  —  AiU^  —  BiUjj  —  AzU^^  —  BzU^jj  —  ^4$  +  Bz^  U^,j 

=  At[-Utt-\■{P  +  B^^-BBl-B^B)u,r, 

+  (4iA  +  BiAi  -  AiB  -  BAi  -  ABi  -  BiA  +  AB  +  Ba) 

+  {piBz  +  BzB\  —  BBz  —  BzB^j  U^nr)  + 

-f  {^BzB\  -}-  BiBz  ~  Bz-B  —  BBz  +  BiAz  -b  AzBi 
•  — BAz  —  AzB  .S24i  *b  A1B2  ~  BzA  —  AB^ 

-{-  {BzBz  +  BzB^  +  BzAz  +  43,82) 

-b  -b  43  -b  4282  -b  B2A2  +  Az-Bz  +  BzA^ 

-b  ^434i  -b  4i43  —  434  —  443  -b  4183  -b  BzA\ 

— ABz  ~  BzA  -b  4381  -b  81 42  ~  AzB  —  BA^ 

-b  ^4343  -b  4343  -b  AzBz  -b  BzAz^ 

-b  (4"  -b  4?  -  44,  -  4i4)  8« 

-b  ^4i42  -b  424i  —  442  ~  42.4)  -b  428(f^^^J  (5.36) 

Examination  of  Eq  5.9  reveals 

Uu  =  (82-b8?-M, -8,8)t/„, 

-b  [AiBi  -b  8,4,  -  4,8  -  84,  -  48,  -BiA  +  AB  -b  8/1) 

+  {pxBz  +  828,  —BBz  —  828)  -b  B^Ujjjjrj,] 
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+  {BzB\  +  B1B3  —  B3B  —  BB:i  4-  BiA'i  +  AzB\ 

—BAz  ~  AzB  +  A1B2  —  B2A  —  AB^ 

+  (^2-63  A-  BzB2  A-  B2AZ  +  AzB-^  ^inm 
+  ^jB|  +  i43  +  A2B2  +  B2A2  +  Az-Bz  +  BzA^ 

+  ^^3^1  +  AiAz  —  AzA  —  j4j43  +  A\Bz  +  BzAi 
—ABz  ~  BzA  +  A2B1  +  B1A2  ~  A2B  —  BA-^ 

+  i^AzAz  +  AzAz  +  AzBz  +  BzA^ 

+  [a^  a  AI-  AAi  -  AiA)  % 

+  ^AiAz  4"  AzAi  —  AA2  —  A2A  )  ^C€«  +  ^2^«€«  (5.37) 

Upon  substitution  of  Eq  5.37  into  Eq  5.36,  the  right-hand  side  of  Eq  5- 36  be¬ 
comes  0  [At^,  A(^^,  At;^]  and  a  second-order  accurate  algorithm  is  assured.  Thus, 
the  scheme  given  by  Eqs  5.21-5.23,  5.33,  and  5.34  is  a  representation  of 

Ut  -f  -f  IJxBrj  +  -j-  TlyGr^  ~  ^xBy^  +  TJ^F "b  ^yGy^  -f  VyGvt,  . 

AO[Ae,Ae,^V^] 

This  scheme  is  hereafter  referred  to  as  ATNSC2,  which  is  shorthand  for  2nd-0rder 
AFIT  TVD  Navier-Stokes  Code.  A  description  of  the  ATNSC2  routines  is  contained 
in  Appendi.x  C.  Also  included  is  a  relative  performance  evaluation  of  the  routines, 
obtained  using  the  CRAY  FLOWTRACE  option. 


5.S.3  Beam-Wai'ming  Algorithm.  The  implicit  Beam- Warming  algorithm  as 
implemented  by  Visbal  [47]  is  used  herein  as  a  comparison  against  the  solutions 
provided  by  .ATNSC.  The  scheme  solves  the  strong  conservation  form  of  the  Navier- 
Stokes  equations.  The  scheme  is  written  in  delta  form  as 


At 

'1+O2 


f  r  ,  O.At  \dA  f  T  ,  g.At  f3/}  atiA'  - 

+  1+O2  la€  de  ]  i  1+O2  la.}  dn-^l  ^  - 

(A  {F  -  A)”  +  I;  (G  -  G.)"]  +  fit  {ftA"-'F„  +  ^A"-'G„,) 


(5.39) 
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where 


A"a  = 

and 

i„  =  dhldu^ 

Fy^  and  G„,  are  evaluated  from 

F4U,U^,U,)  =  F 
G,{U,U^,U,)  =  G 

For  steady-state  calculations,  first-order  accurate  Euler-implicit  time  differencing  is 
utilized  by  setting  0i  =  I  and  Oo  =  0  .  Second-order  time  accuracy  is  achieved  by 
setting  =  1/2  and  02  =  0  to  obtain  a  trapezoidal  time  differencing.  Trapezoidal 
differencing  is  used  for  all  computations  where  the  a  time-dependent  solution  is  of 
interest.  Eq  5.39  is  implemented  using  second-order  appro.vimations  for  the  spatial 
derivatives. 

To  maintain  numerical  stability  and  provide  smooth  solutions,  e.xplicit  fourth- 
order  damping  is  added  to  the  right-hand  side  of  Eq  5.39  as 


:  6^”+!  -  6'” 


(5.40) 


Gv  =  {VtF^  +  V’jGv)  I J 

By  =  dGy/dUy 


(5.41) 


{UM^)-^Fy,{U,Ur,) 

'v.  {UM  +  Gy,{UM 


(5.42) 


-c4AtJrl5lUr:. 


(5.43) 


and 


-^;Au-j6iurj 


(5,4.1) 


Implicit  second-order  damping  is  inserted  with  respect  to  the  implicit  operators  as 


-u^AUr;SlJi.jI  (5.45) 

and 

-w^AUrJS-A,l  (5.46) 
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In  Eqs  5.43  through  5.46,  is  a  second-order  accurate,  central-difference  operator 
used  to  approximate  the  derivative  with  respect  to  /.  The  nominal  values  of  the 
damping  coefficients  are  [27] 


wf  =  0.02  u’J  =  0.04 
w?  =  0.25  ujf  =  0.25 


(5.47) 
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VI.  Viscous  Results  and  Conclusions 

6.1  Burgers’ Equation 

To  provide  a  performance  comparison,  the  Lax-VVendroff  and  ATNSC  schemes 
are  applied  to  the  linearized  version  of  the  viscous  Burgers’  equation 

Ut  +  CUx  =  flUxx  (6.1) 

with  periodic  boundary  conditions 

u{0,t)  =  n{27r,l)  (6.2) 

and  the  initial  condition 

n(.T,0)  =  csin(A-x’)  (6.3) 

Eq  6.1  has  the  time-dependent  solution 

U {x,  t)  =  sin{/:(.x  -  cl)]  {6.-i ) 

The  equations  can  be  non-dimei).«ionalized  using 

X  =  kx- 

u  =  u/c  (6.5) 

L  =  Ike 

to  obtain 

n(x,0)  =  sinx  (6.6) 

LLi-lV)  =  (x  —  <) 
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where 

Re  =  cKf^ik)  (6.7) 

Rather  than  perform  a  stability  analysis  for  the  ATNSC  scheme  with  the  viscous 
term  added,  the  exact  stability  condition  for  the  La.x-VVendrolF  method  is  u.sed  [l]. 
In  non-dimensional  form  this  limitation  becomes 


F'igures  6.1  and  6.2  show  the  results  obtained  using  first-order  time  accurate, 
'P”  =  (1/Re)«^  ,  versions  of  the  La.x-VVendroff  and  ATNSC  schemes.  This  is  similar 
to  the  method  used  in  the  implicit  application  of  TVD  algorithms  to  viscou.s  ecpia- 
tions  and  is  utilized  herein  to  demonstrate  the  necessity  for  second-order  accuracy 
at  low  Reynolds  numbers. 

Solutions  are  initially  computed  under  a  CFL  restriction  of  0.80  with  Re  = 
10000  using  49  (Case  I),  99  (Case  II),  and  199  (Case  III)  cells  respectively.  The 
solution  was  carried  out  to  i  =  56.665  .  Figure  6.3  shows  the  same  schemfjs  under 
a  CFL  restriction  of  0.95  and  Re  =  10000  using  49  cells  only.  The  norm  given  witli 
each  plot  is  defined  as 

IIm  -  U.\\  =  -  £1'^  (fi-9) 

H  is  the  exact  solution  of  Eq  6.6  shown  as  the  .solid  cui  ve  in  the  figures. 

Figure  6.1  shows  that  the  Lax-VVeudrolf  scheme  exhibits  a  phase  shift  fur  1, 
which  diminishes  as  the  number  of  cells  is  increased.  Given  the  smooth  initial  data  of 
Eq  6.6,  the  overall  performance  of  the  Lax- Wendroff  scheme  is  lathcr  good.  However, 
the  phase  shift  a5.sociated  with  the  Lax- Wendroff  scheme  does  not  appeal  with  the 
ATNSC  .schemes. 

Figure  6.2  shows  the  results  obtained  u.sing  the  first-order  .ATN.SC  scheme.  The 
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majority  of  the  error  is  concentrated  in  the  areas  of  rapidly  changing  gradients.  Note, 
however,  that  the  ATNSC  norms  are  smaller  in  every  case  than  the  Lax-Wendroff 
norms. 

Figure  6.3  shows  the  effect  of  increasing  the  CFL  restriction  to  0.95  for  the 
two  schemes.  The  increase  in  the  CFL  (or  Courant)  number  tends  to  lessen  the 
phase  shift  associated  with  the  Lax-VVendroff  scheme,  while  allowing  the  computed 
solutions  of  ATNSC  to  more  closely  match  the  e.xact  solution  in  the  regions  of  rapidly 
changing  gradients. 

Figure  6.4  displays  the  results  of  a  second-order  version  of  the  .ATNSC  scheme 
applied  to  Eq  6.6  at  Re  =  10000.  Second-order  accuracy  is  arrived  at  through  the 
use  of 

K  (6-10) 

which  is  consistent  with  Eq  5.31.  A  CFL  restriction  of  I  was  used  with  no  observed 
difficulty  and  the  norms  were  between  one  and  two  orders  of  magnitude  lower  than 
either  of  the  first-order  schemes. 

Results  of  a  more  severe  test  of  the  first-order  version  of  .ATNSC  are  shown  in 
Figure  6.5.  The  first-order  algorithm  was  applied  K.*  the  viscous  Burgers'  equation 
with  Re  =  100  and  i  =  25.  The  lower  Reynolds  number  represents  an  increase  in 
significance  of  viscous  forces  over  the  previous  case.  The  figure  clearly  shows  the 
degradation  of  the  solution  near  the  extrema,  as  well  as  small  oscillations  present  in 
the  solution.  The  CFL  number  had  to  be  reduced  as  spacing  was  reduced  in  order 
to  maintain  stability.  In  fact,  the  CEl-  numbers  used  were  apparently  on  the  edge 
of  the  stability  limit. 

Solutions  using  the  second-order  accurate  version  of  .AT.NSC  for  Re  =  100  are 
shown  in  Figure  6.6.  The  solution  is  crisp  near  the  peaks  and  is  free  of  oscillations. 
The  norms  are  an  order-of-magnitude  lower  than  those  of  the  first-ordei-  scheme. 
Finally,  the  .second-order  scheme  wa.s  utilized  with  a  CFL  number  of  i  for  each  cell 
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size.  The  resuits,  shown  in  Figure  6.7,  are  indistinguishable  from  those  of  Figure  6.6 
except  that  the  increase  in  CFL  number  further  decrea.sed  the  norms. 

Overall,  all  the  first-order  schemes  perform  rather  well  for  the  viscous  Burgers’ 
equation  with  periodic  boundary  conditions,  smooth  initial  data,  and  high  Reynolds 
number.  However,  the  Lax-VYendroff  solution  can  be  expected  to  show  the  same  type 
of  oscillations  as  for  Riemann’s  problem  if  smooth  initial  data  is  not  specified.  It 
should  also  be  noted  that  the  first-order  TVD  scheme  performs  best  when  the  CFL 
restriction  is  as  close  to  1  as  possible. 

The  situation  is  not  so  favorable  for  the  first-order  TVD  scheme  at  lower 
Reynolds  numlier.  The  first-order  scheme  exhibits  decreased  accuracy  and  severe 
stability  restrictions.  Thus  the  first-order  scheme  is  wholly  unsuited  to  the  calcula¬ 
tion  of  unsteady  flows  at  low  Reynolds  numbers. 

The  behavior  of  the  second-order  TVD  scheme  is  very  encouraging.  Accuracy 
is  superb  and  the  scheme  is  very  robust  where  stability  is  concerned.  CFL  numbers  as 
high  as  l.l,  in  relation  to  Eq  6.8,  were  tested  with  no  instability.  The  Lax-VVendroff 
time  step  limifation  given  by  Eq  6.8  is  ob\iously  somewhat  over  restrictive  in  this 
ca.se. 
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Figure  6.2.  ls(.-Orclcr  .'Vl’N.SC  Scheme  .Applied  to  Viscous  Burgers'  Fcpiation. 
CFL=0.S 
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Figure  6.3.  Ist-Order  Lax-VVendrofr  and  .VP.X.SC  .Scheme.-?  .Applied  to  V’iscout.  Bnrg- 
er.s'  Equation,  CFL=0.95 
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Figure  6.4.  2nd-Orcler  .AT:\'.SC  .Scheme  .Applied  to  Vi.scou.s  Burger 
CFL=1.0 


6.2  Bouiidanj  Conditions  for  Viscous  flow 


Boundary  conditions  for  viscous  flows  are,  in  general,  more  straight  forward 
than  their  inviscid  counterparts  described  in  Section  3.2.  .At  the  wall,  the  inviscid 
surface  tangency  condition,  Eqs  3.11  and  3.10,  is  replace  by  the  viscous  no-slip 
requirement: 

«  =  0 

(6.11) 

e  =  0 

Simplified  wall  temperature  contition.s  representing  either  an  adiabatic  wall 

r...  =0  (6.12) 


or  constant  temperature  wall 


(6.13) 


arc  used  in  the  current  study,  depending  on  the  flow  of  interest.  With  the  wall 
mapped  to  a  constant  77  coordinate,  the  pressure  at  the  wall  is  obtained  by  solving 
the  normal-momentum  equation: 


Pn  yjnl  +  p'l  =  {fxpr  +  )  Vi.  +  ('/j  +  '/y)  Prt 

+'/!,  (2/t  d-  A)^.  -r  iPj  (2/f  -r  A),J  } 

d-  (4«4  d-  7/y»„)  {px  iiyl'i  d-  ilyP,,)  -b  Vy  +  Vxl'v)] 

d-/«/r  d-  2<,//,,»^.,  -r 

-}-(2/t  -f  A);/;,  2^y/7.;7v,j  Hh  77'7-,j,,^ 


Flow  at  the  inlet  and  c.vit  of  the  compul.alional  domain  is  a.ssumcd  to  be 
inviscid.  Inflow  and  outflow  relations  from  Section  3.2  are  thus  used  to  determine 
flow  quantitic.s  at  thc.se  boundaries.  .\s  staled  in  .Section  3.2,  for  supersonic  outflovv 
all  quantities  must  be  e.x'trapolated  from  the  interior  of  the  domain.  In  prad  ce. 
this  extrapolation  is  also  performed  in  (he  .Mib.sonic  boundaiy-layei  emfjedded  in  the 


supersonic  outflow.  For  the  cases  to  be  considered  herein,  no  adverse  eflects  of  this 
e.xtrapolation  are  noied. 

6.3  Shock- Boxindary  Layer  Interactioii 

An  indepth  e.xperimcnt  in  laminar  shock-boundary-layer  interaction  was  car¬ 
ried  out  by  Hakkinen  et  al.  [IS]  in  1959  at  the  Massachusetts  Institute  of  Technology 
under  the  sponsorship  of  the  National  .Advisory  Committee  for  .Aeronautics.  De¬ 
tailed  measurements  were  made  of  pressure  distribution,  skin  friction  coefficient, 
and  velocity  profiles  for  a  number  of  combinations  of  overall  pressure  ratio,  p/fpcc-, 
and  shock  Reynolds  number,  Rcx,,  at  a  freestream  .Mach  number  of  2.0  for  a  shock 
wav’e  impinging  upon  a  flat  plate  boundaiy- layer.  The  most  recognizable  of  these 
in  the  CFD  community  is  the  case  of  Figure  6b  of  reference  [IS],  The  overall  pres¬ 
sure  ratio  for  this  case  is  1.40  at  a  shock  Reynolds  number  of  2.96  x  10®,  based  on 
a-j  =  4.978  cm.  It  was  pointed  out  by  Dcgrez.  Boccadoro,  and  Wendt  [llj  that  this 
has  been  used  as  a  test  case  by  numerous  researchers  (Skoglund  and  Gray  in  1969; 
MacCormack  in  1971  and  1982;  Hanin,  Wolfshtein,  and  Landau  in  1974;  Beam  and 
Warming  in  1978;  and  Dawes  (10]  in  1983).  Lion  also  used  this  as  a  test  case  as 
recently  as  1989  f26|. 

The  e.xperimental  pressure  and  skin  friction  profiles  for  this  case  arc  shown  in 
Fiqure  6.8.  and  a  sketch  of  the  wave  structure  is  shown  in  Figure  6.9.  The  friction 
coefficient.  Cj.  is  defined  as 


where  r^.  is  the  normal  component  of  shear  stress  at  the  wall 

=  (6.16) 
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and  f/eo  is  the  dynamic  pressure 


<ZcO  — 

With  the  tangential  velocity  given  by  Eq  3.11,  and  the  wall  mapped  to  a?/  =  constant 
coordinate,  can  be  written  as 


—  P  {Vy^71  Vx^T]')  (6. lb) 

No  negative  values  of  skin  friction  are  shown  because  the  total-head  tube  was  not 
able  to  reliably  indicate  negative  shear  values.  Locations  where  the  e.xperimental 
skin  friction  may  have  been  negative  are  shown  by  downward  pointing  arrows  in 
Figure  6.8.  While  the  accuracy  with  which  the  pressure  profile  can  be  calculated  has 
greatly  improved  since  MacCormack’s  calculations  [28],  there  has  been  essentially  no 
progress  in  matching  the  skin  friction  profile.  This  includes  the  overall  shape  of  the 
profile  as  well  as  the  location  of  the  separation  and  reattachment  point.s.  MacCor¬ 
mack’s  calculations  failed  to  show  the  characteristic  plateau  in  the  pressure  profile, 
and,  while  obtaining  a  fair  prediction  of  the  separation  point,  he  predicted  reat- 
lachment  ahead  of  the  experimental  data.  In  addition,  the  friction  coefficient  after 
reattachment  is  approximately  20%  lower  thair  that  suggested  by  the  experimental 
data.  Liou  [26]  obtained  a  fair  matching  to  the  pressure  profile  but  failed  at  pre¬ 
dicting  the  skin  friction  profile  in  the  regions  of  adverse  pressure  gradient.  In  eveij 
case  known  to  the  author,  eveir  those  that  somehow  managed  to  accurateh  ]rredict 
separation  and  reattachment  points,  the  ultimate  skin  friction  lc\el  after  reattach¬ 
ment  remains  18%  -  20%  low.  Liou  goes  .so  far  as  to  state  that  this  discrepancy  in 
the  skin  friction  level  may  be  due  to  transition  of  the  boundary-layer  from  laminar 
to  turbulent  immediately  in  the  interactior  region  (26).  This  is  in  direct  contrast  to 
the  experimental  velocity  profiles  of  reference  [18]  and  contrary  to  the  observations 
of  the  experimenters  [18]. 
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Figure  6.8.  Experimental  Pre.ssure  and  Skin  Friction  Profiles 


INCIDENT 


Figure  6.9.  Flosvfiold  Structure 


Figure  6.10.  Grid  Used  in  Shock-Boundary  Layer  Interaction  Investigations 


The  grid  used  for  the  numerical  investigations  is  shown  in  Figure  6.10,  with 
133  points  in  the  axial  direction  and  60  points  in  the  normal  direction.  Spacing  is 
held  constant  in  the  axial  direction  at  Ax/xshock  =  0.013  and  ranges  in  the  normal 
direction  from  an  initial- value  at  the  wall  of  ^tj/xshock  =  6.78  x  10“'’  to  a  final  value 
olAij/Xshock  =  I -12  X  10"'^  at  the  upper  edge.  Grid  densities  are  chosen  comparable 
to  those  used  l)y  MacCormack  (28),  Dawes  [10],  and  Lion  (26)  to  provide  a  comparison 
based  on  similar  grids. 

6.3.1  .A'TNSC  Solutions.  The  computational  domain  is  initialized  at  the  uni¬ 
form  freestream  conditions  to  the  left  of  the  point  along  the  upper  boundary  at 
which  the  shock  is  generated.  Post-shock  conditions  aie  applied  downstream  of  this 
point.  .An  adiabatic  wall  condition  is  used  to  obtain  the  wall  temperature  along  the 
plate  and  the  nomal  momentum  equation  is  solved  to  obtain  the  wall  pre.ssure  in 
combination  with  the  no-slip  velocity  constraint  at  the  wall. 

Figures  6.11-6.22  show  the  results  of  applying  the  .ATN.SC  algorithms  to  this 
test  case.  The  data  repre.sented  by  the  figures  can  be  taken  to  be  the  solution  pro- 
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vided  by  either  ATNSCl  or  ATNSC2,  since  at  this  Reynolds  number  both  algorithms 
provided  exactly  the  same  results. 

Figures  6.11  and  6.12  depict  the  solution  obtained  with  e  =  0  in  the  nonlinear 
fields  and  e  =  0.05  in  the  linearly  degenerate  fields.  Note  that  this  is  in  contrast 
to  the  values  of  e  =  0(0.1)  for  the  nonlinear  fields  and  e  =  0  for  the  linearly  de¬ 
generate  fields  typically  utilized  for  the  inviscid  calculations  of  Part  I.  Values  of  t 
up  to  0.025  for  the  nonlinear  fields  were  found  to  have  no  noticeable  effect  on  the 
solution  while  the  e  value  used  in  the  linearly  degenerate  fields  significantly  alter 
the  solution,  as  will  shortly  become  apparent.  The  pressure  profile  of  Figure  6.11 
clearly  shows  the  pressure  ri.se  to  scpaiation,  the  constant  pressure  plateau  within 
the  separated  region,  and  the  pressure  rise  to  leattachment  as  described  in  refer¬ 
ence  (18).  The  most  noticeable  aspecl  of  the  pressure  profile  is  the  slightly  lower 
value,  as  compared  to  the  experimental  data,  within  the  separation  region.  The 
reason  for  this  is  unknown,  although  the  trend  was  consistent  throughout  the  invesi- 
gation.  The  skin  friction  profile  of  Figure  6.U  contains  several  regions  of  interest. 
-First,  there  is  a  very  slight  oscillation  in  the  friction  coefficient  leading  up  to  the 
sharp  drop  just  prior  to  separation.  This  was  observed  for  values  of  >  0.025, 
and  in  fact,  skin  friction  was  severly  oscillatory  at  c>.,i  =  0(0.1)  which  are  not  un¬ 
usual  values  when  £2.-1  7^  0  for  inviscid  calculations.  The  length  of  the  separation 
region  was  underpredicted  in  that  delay  ed  .sepaiation  and  premature  leatlaclnnent 
were  observed.  This  again  appears  to  be  an  artifact  of  the  values  of  t2..i  used,  as 
will  become  apparent  upon  examination  of  .•'ub.seciuent  figures.  The  skin  friction 
profile  beyond  reattachment  shows  a  rapid  ii.se  to  the  ultimate  value,  although  the 
ultimate  value  shows  much  better  agreement  with  the  experimental  data  than  that 
obtained  through  previous  investigations  known  to  the  author.  Figure  6.12  provides 
a  visualization  of  the  wave  structure  through  50  equally  spaced  pressure  contours 
between  the  upstream  and  downstream  prcssuics.  The  AT.N’SC  algorithm  provides 
high-resolution  capturing  of  all  the  pertinent  flow  structures.  The.se  include  the  gen- 


crated  shock,  the  leading-edge  shock  and  accompanying  expasion,  separation  shock, 
expa.nsion  fan,  and  reattachment  shock. 

The  values  of  e2.‘i  were  lowered  to  0.02.5  for  the  solution  depicted  in  Fig¬ 
ures  6.13  and  6.14.  The  pressure  profile  of  Figure  6.13  is  identical  to  that  of  Fig¬ 
ure  6.11  except  for  the  extension  of  the  constant  pressure  plateau  slightly  upstream 
and  downstream.  This  is  due  to  the  inciease  in  the  length  of  the  separation  region 
apparent  upon  comparison  of  Figures  6.11  and  6.13.  Calculated  separation  and  reat¬ 
tachment  points  agree  extremely  well  with  the  experimental  data.  Note  also  that 
there  is  no  longer  aii\  oscillation  in  the  friction  coefficient  in  the  upstream  region  and 
that  the  rise  to  the  ultimate  downstt earn  friction  value  is  more  gradual.  This  is  the 
first  numerical  solution  known  to  the  author  that  correctly  predicts  the  separation 
and  reattachment  points  as  well  as  the  correct  downstream  friction  coefficient.  Com- 
.pu ted  pressure  contours  for  this  particular  case  are  presented  in  Figure  6.14.  Wave 
structure  is  very  similar  to  that  of  Figure  6.12  except  for  the  enhanced  structure  in 
the  interaction  region,  due  to  the  lengthening  of  the  separation  region.  .-\  lengthed 
separation  region  also  provides  enhanced  re.solution  of  the  expansion  fan  in  that  it 
is  not  so  tightly  packed  between  the  shocks. 


Values  of  c  necessary  to  produce  an  acceptable  solution  are,  as  alluded  to 
previously,  an  order  of  magnitude  smaller  than  the  values  that  arc  often  used  for 
inviscid  flow.  Since  the  vast  majority  of  viscous  T\'D  research  has  been  conducted 
for  hypersonic  flows,  an  answer  was  sought  in  the  appropriate  literature.  Exami¬ 
nation  of  references  (24],(32|,  and  (34)  revealed  that  values  of  0.05  <  c  <  0.25  were 
commonly  used  for  hypersonic  flows  in  the  range  1  <  <  25  with  c  as  high  a.s 

0.5  in  some  instances.  However,  it  wa.s  discoxcred  that  these  re.seaichers  were  using 


variable  i.sotrcp'C  damping  attributed  to  Vee  [13]  and  anisotropic  damping  due  to 
Martinelli  (29).  In  the  normal  direction,  isotropic  damping  is  applied  to  the  nonlinear 
fields  as 


e"  =  cAt 


\u  ■  V^(  -r  )u  •  V//)  4-  -  (/j^  -f  k„) 


(6.19) 
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and  in  the  axial  direction,  anisotropic  damping  is  applied  to  all  fields  as 


(6.20) 


where  =  |«  •  V)t|  +  c|V^•|. 

These  changes  were  incorporated  into  a  version  of  the  ATNSC  code  and  a 
solution  calculated  again  using  c  =  0.0  in  the  nonlinear  fields  and  e  =  0.0-5  in  the 
linearly  degenerate  fields.  This  particular  solution  is  depicted  in  Figures  6.15  and 
6.16.  Figure  6.15  shows  that  this  solution  is  identical  to  that  of  Figure  6.11  except 
for  a  small  'hange  in  the  skin  friction  profile  near  the  minimum  value.  Examination 
of  Figure  6.16  reveals  a  sharper  resolution  of  the  wave  structure  than  Figure  6.12, 
more  in  line  with  the  structure  of  Figure  6.14. 

Lowering  the  value  of  c  in  the  linearly  degenerate  fields  to  0.025  resulted  in  the 
solution  of  Figures  6.17  and  6. IS.  .At  this  lower  value  of  c,  the  profiles  and  contours 
of  Figures  6.17  and  6.18  appear  identical  to  those  of  the  constant  damping  case, 
Figures  6.1-3  and  6.14. 

Based  on  this  set  of  solutions,  it  appears  that  the  change  in  skin  friction  pro¬ 
files  and  wave  structure  with  c  is  not  a  strong  I'uncticn  of  variable  versus  const.uU 
damping.  Research  was  then  conducted  into  whether  anyone  had  observed  the  same 
phenomenon  in  a  similar  flow  regime.  Seider  and  Ilanel  (39)  have  recentlv  observed 
similar  phenomenon  in  regards  to  traiusoiiic  airfoil  drag  prediction.  Tliev  bimulatod 
the  transonic  flow  about  a  R.*\E  2822  airfoil  {M^  =  0.73. a  =  2.79°)  using  several 
TVD  schemes  applied  to  the  thin  laver  N'a\  icr-Slokes  equations,  including  a  .'scheme 
based  on  Roe's  cipproximate  Riemann  .■solver.  They  found  that  the  Roe  ba.scd  scheme 
provided  the  best  overall  results,  but  that  a  certain  sensitivity  to  the  values  of  c  for 
the  linearly  degenerate  fields  existed  for  all  their  schemes.  They  used  values  of  t  —  0. 
0.1,  and  0.2  and  found  that  c  =  0.2  resulted  in  a  4%  decrea.se  in  drag,  while  leaving 
lift  unchanged.  This  is  consistent  with  the  changes  in  skin  friction  while  pre.ssure 
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Figure  6.17.  Pressure  and  Skin  Friction  Profiles  (ci  =  £3  =  O.e^ 


Figure  6. IS.  Pressure  Contours  (cj 


remains  constant  <as  seen  herein.  They  then  examined  a  flat  plate  at  M^o  =  0.5  and 
Ret  =  5000  and  found:  that  shin  friction  in  this  case  increased  with  increasing  t. 
Finally,  they  doubled  the  number  of  grid  points  in  the  boundary-layer,  from  7  to  14, 
and  found  that  this  totally  removed  the  t  dependence.  It  appears  that  this  behavior 
is  common  to  flows  in  the  transonic  and  low  supersonic  regimes  and  the  effect  of  t 
must  be  analyzed  whenever  a  solution  is  computed  at  these  Mach  numbers. 

Three  final  solutions  using  .4TNSC  are  presented  showing  the  above  mentioned 
behavior.  First,  the  variable  damping  algorithm  is  used  to  arrive  at  the  solution  of 
Figures  6.19  and  6.20  using  c  =  0.025  for  all  fields.  This  solution  is  identical  to  that  of 
Figures  6.13  and  6.14,  thus  supporting  the  a,ssertion  that  the  value  of  c  in  the  linearly 
degenerate  fields  is  the  primary  influence.  Halving  t  led  to  the  computed  solution 
shown  in  Figures  6.21  and  6.22.  .-\gain,  the  pressure  profile  remains  essentially  the 
same  as  all  other  cases,  but  the  skin  fiiclion  levels  have  decreased  slightly,  resulting 
in  premature  separation  and  delai^ed  reattachment.  Finally,  the  number  of  grid 
points  in  the  boundary-layer  was  doubled,  from  10  to  20.  Solutions  are  presented  in 
Figure  6.23  for  c  values  of  0.0125.  0.025,  and  0.035  using  this  new  grid.  The  pressure 
profile  remains  unchanged  except  for  a.  decrease  in  the  length  of  the  pressure  plateau. 
Skin  friction  changes  only  slightly  upstieani  of  the  interaction  region,  but  dro|)s  to 
zero  more  rapidly  than  is  the  case  in  Figure  G.  13.  Separation  and  reattachment  points 
are  correctly  predicted,  and  the  rise  to  the  final  skin  friction  level  more  closely  follows 
the  experimental  data  than  that  of  the  previous  solutions.  The  e  dependence  has 
been  removed,  within  the  range  0.0125  <  <  <  0.035,  consistent  with  the  ob.servartions 
of  Seider  and  Hanel  [39]. 

f 

G.3.2  Beam-Wanning  Sohil.ion.'>.  The  approximate-factorization  algorithm  of 
Beam  and  Warming,  Eq  5.39  as  implemented  by  Visbal  (47),  is  applied  to  this  test 
case  as  a  comparison  against  the  .AT.\'SC  algorithm.  Figures  6.24  and  6.25  de|)ici  the 
Beam- Wanning  solutions  using  the  nominal  recommended  value  of  the  second  and 
fourth-order  damping  coefficients  [27].  A  value  of  0.25  i.s  used  for  the  second-ouier 
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Figure  6.21.  Pressure  and  Skin  Friction  Profiles  (ci  =  ej  =  =  0.0125) 


Figure  6.22.  Pressure  Contours  (ei  =  62  =  £3  =  e-i  =  0.0125) 
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Figure  6.23.  Pressure  and  Skin  Friction  Profiles  (20  Points  in  the  Boundary-Layer) 

coefficients  in  both  the  axial  and  normal  directions  while  a  value  of  0.02  is  used  for 
the  fourth-order  coefficient  in  the  axial  direction  and  0.04  for  the  normal  direction. 
The  pressure  profile  of  Figure  6.21  is  very  similar  to  the  ATNSC  profiles  except  for 
a  slightly  shortened  length  of  the  pressure  plateau  in  the  separation  region.  The 
skin  friction  profile  is  similar  to  the  profiles  appearing  in  the  literature  for  various 
algorithms  applied  to  this  problem.  Beam-Warming  predicts  early  separation  while 
the  reattachment  point  is  in  good  agreement  with  the  experimental  data.  Most 
noticeable  is  the  under-prediction  of  the  skin  friction  level  bevond  icattacluneiit . 
approximately  1.5%  below'  the  experimental  value.  Figure  6.25  depicts  a  significant 
decrease  in  re.solution  of  the  w'ave  structure  as  compared  to  the  ATNSC  solutions. 
Shock  waves,  compression  waves,  and  expansion  waves  are  all  smeared  to  a  much 
greater  extent  than  those  obtained  with  TVD. 

Figures  6.26  and  6.27  show  the  .solution  computed  when  the  damping  coeffi¬ 
cients  were  double  tho.se  of  the  pro\  ions  solution.  The  pressure  profile  shows  a  fui  tlici 
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Figure  6.24.  Pressure  and  Skin  Friction  Profiles 
(w«  =  0.02,  w?  =  0.04,  w?  =  wf  =  0.25) 


Figure  6.25.  Pre.ssure  Contours  (lc^  =  0.02, =  0.04, 
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reduction  in  the  length  of  the  plateau  within  the  separation  region.  Predicted  sepa¬ 
ration  and  reattachment  points  are  in  good  agreement  with  the  experimental  data, 
while  the  ultimate  friction  level  is  still  below  that  of  experiment.  There  is  no  oscil¬ 
lation  of  the  skin  friction  profile  in  this  highly  damped  case.  The  pressure  contours 
of  Figure  6.27  are  more  highly  smeared  than  those  of  Figure  6.25,  but  the  increased 
damping  results  in  a  reduction  of  the  oscillatory  effects  upstream  and  downstream 
of  the  shocks,  compressions,  and  expansions. 

Finally,  Figures  6.28  and  6.29  depict  a  solution  using  half  the  damping  values 
of  Figures  6.24  and  6.25.  The  profiles  of  Figure  6.28  show  that  the  pressure  plateau 
in  the  separation  region  has  lengthened  over  that  of  Figure  6.26,  more  in  line  w'ith 
Figure  6.24,  but  stiil  is  less  in  extent  that  that  of  the  ATNSC  solutions.  .Again, 
early  separation  is  evident  but  reattachment  is  in  line  udth  the  experimental  data. 
Oscillations  in  the  friction  profile  downstream  of  reattachment  have  reappeared,  and 
the  ultimate  skin  friction  remains  approximately  15%  below  the  experimental  value. 
Wave  structure,  Figure  6.29,  is  much  less  smeared  than  in  Figure  6.27,  but  shows 
the  same  oscillatory  effects  as  Figure  6.25. 

The  solutions  obained  using  the  .ATNSC  TVD  algorithm  clearly  demonstrate 
that  it  is  possible  to  obtain  very  accurate  estimations  of  separation  and  reattachment 
points  while  at  the  same  time  maintaining  high  degrees  of  resolution  of  the  u.ive 
structure.  The  Beam-Warming  solutions  are  consistent  with  those  in  the  liteiatuie 
for  numerous  algorithms.  Only  the  T\’D  based  .ATNSC  algorithm  incorporates  the 
necessary  physical  constraints  to  achieve  high  accuracy  and  resolution. 

.ATNSC  .solutions  are  obtained  using  a  constant  CFL  number  of  0.95,  under  the 
timestep  restriction  of  Eq  1.28.  The  solution  is  monitored  until  no  change  is  olxsetwed 
in  the  skin  friction  profile,  typically  requiring  4000  time  steps  to  achieve  stead} - 
state  convergence.  The  data  processing  rate  for  ATNSCl  with  constant  damping  is 
1.6071  X  10"’  seconds  per  grid  point  per  time  level  for  the  CRAY  X-.MP/216.  Data 
processing  rates  for  the  other  .AT.NSC  versions  arc  contained  in  .Appendix  C.  Beam- 
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Figure  6.26.  Pressure  and  Skin  Friction  Profiles 
(w<  =  0.04,  =  0.08,  =  a;?  =  0.50) 


Figure  6.27.  Pressure  Contours  (a;|  =  0.04,0; 
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Warming  solutions  are  obtained  using  a  CFL  number  of  5.0,  again  monitoring  until  no 
change  is  observed  in  the  skin  friction  profile.  This  typically  requires  1000  iterations 
for  Beam-Warming  to  achieve  steady-state  convergence.  The  Beam- Wanning  data 
processing  rate  is  1.9316  x  10"®  seconds  per  grid  point  per  iteration  for  the  CRAY 
X-MP/216. 


6-4  Unsteady  Shock-Indxiced  Heat  Transfer 

6.4.1  Heal  Transfer  Theory  of  Mirels.  The  numerical  solutions  to  the  un¬ 
steady,  shock-induced  heat  transfer  problem  are  compared  herein  against  the  lheoi\ 
of  Mirels  (30,  31]  as  implemented  by  .Schlichting  [37].  The  heat  flu.x  at  the  wall  can 
be  written  as 

<l  =  hATa,u-r,,)  (6.21) 

where  the  adiabatic  wall  temperature,  Taw.  i-?  given  by 


(1  + 


(6.22) 


and  the  local  convective  heat  transfer  coefficient,  /(,  is  defined  as 


,  :\'u,h 


(6.23) 


For  the  laminar  flow  of  a  perfect  gas  u  ith  weak  shocks,  Mirels  pun  ides  the  following 
appro.ximation  for  the  recovery  factor,  r 


=  Pr" 


(6.2-1) 


where 


o  =  O..39  - 


0.05 


1 


(6.25) 
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and  Us  is  the  shock  propagation  speed.  The  local  Nusselt  number  tor  an  unsteady 
flow  is  defined  as 

Nus..  =  0.oCfRe  Pr^  (6.26) 


.vhere  Cf  is  the  local  skin  friction  coefficienl 


Cf  =  ( 1  -  0.346'^’"' 


[!  V  n 


VluV 


U. 


(6.27) 


Rt  is  the  Revnolds  number 


Re  = 


(6.28) 


and  thesubsciipt  w  refers  to  propcrtie:  evaluated  at  the  wall  temperature.  Time  in 
Eqs  6.23  and  6.28  is  referenced  to  the  time  of  shock  pa.ssage.  Finally,  the  e.xponent 
of  the  Prandtl  number  in  Eq  6.26  is 


A  =  0.35  + 


0.15 


1  -  U^/Us 


(6.29) 


The  above  relations  are  used  to  calcuh  tlic  theoretical  heat  flux  at  the  plate  wall. 


6.4.S  Nvmerical  Heat  Transfer  Solutions.  The  final  test  case  for  the  .ATNSC 
algorithm  is  the  prediction  of  unsteady,  iaminar  heat  transfei  due  to  a  sliock  wave 
moving  down  a  flat  plate.  The  origin  of  t  hi.s  lest  caase  is  the  work  of  Smith  [41]  who 
used  a  shock  tube  to  study  the  heat  tian.^fei  to  a  sharp-edged  fiat  plate,  crcaiing 
ratios  of  gas  temperature  to  surface  lempeiature  typical  of  those  in  gas  turbine 
engines.  .A  schematic  of  Smith'.s  experimentai  appartus  is  shown  in  Figure  G.30. 

The  case  under  con.$ideralion  here  i.'-  sepresentative  of  data  set  .A  of  Smitli  [4 1  j. 
The  governing  parameters  are  the  shovk  .\lc.ch  number  (Ms)-  pressure  in  the  diiveii 
.section  {P\).  and  temperature  in  the  driven  .sechon  (Ti).  The  wall  tempera* a:.: 
on  the  flat  plate  is  held  constant  in  the  calculations  at  T|.  Using  iV/,  =  1.095, 
P]  =  49102.800  .U/nP.  and  Tj  =  297. 128  f\  re.^ults  in  a  shock  pre.ssure  ratio  of  1.232. 


Figure  6.30.  Flat  Plate  Mounted  in  Shock  Tube 

a  ten.perature  ratio  of  1.062,  a  steady  velocity  behind  the  shock  of  -52.368  iv/s,  a 
steady  ^iach  number  behind  the  shock  of  O.MT,  and  a.  steady  Reynolds  number  be- 
h.ind  the  shock  of  92482.  Shock  Mach  number,  driven  section  pressure,  and  driven 
section  temperature  are  consistent  with  data  set  A  of  reference  [11].  E.vperimen- 
tal  measurements  of  the  shock  Mach  number  are  only  accurate  within  ±2%.  and 
can  significantly  effect  the  level  of  agreement  141]  between  theory,  experiment,  and 
numerical  .solution.  Thu.s,  the  numerical  sulrMons  and  theoretical  values  used  tin, 
nominal  shock  Mach  number  of  1.095  for  comparison.  A  .solution  for  =  1.117. 
2%  above  nominal,  is  presented  as  a  final  compari.son. 

Initial  cond;'  io.is  for  the  computations  consist  of  placing  the  shock  just  ahead  of 
the  plate  at  time  lo  by  establishing  pre-shock  and  post-shock  conditions  on  eithei 
side  of  I  he  point  sclw  ted.  Tiie  shock  is  then  allowed  to  nio\e  freel\  as  time  progre.sses. 
With  the  temperature  held  at  7j  the  normal  .j.eme.itam  equation  is  .solved  to  obtain 
the  pressure  at  the  wall,  in  combination  with  the  no-slip  constiaiul  at  the  wMl.  The 
numerical  solution  is  sampled  at  a  point  5.080  x  10”^/u  downsticam  of  the  leading 
edge  to  obtain  the  heat  flu.x  at  this  point.  This  point  wa.s  chosen  consistent  with 
the  first  sampling  point  of  Smith  in  the  shock  tube  experiment.  The  computations 
t.'L  >-,.i:iied  out  to  a  lime  of  approximately  one  millisecond,  the  api)roximate  time 
ol  transition  to  turbulent  flow  as  noted  by  Smith  [llj.  Two  different  grids  are  u.sed 
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Figure  6.31.  Initial  Grid  for  Hea:  Flu.x  Solutions 

in  this  study,  the  first  of  which  is  shown  in  Figure  6.31.  This  grid  consists  of  201 
points  in  the  axial  direction  and  31  points  in  the  normal  direction.  The  grid  spacing 
is  held  constant  in  the  axial  direction  at  A.r  =  2. .510  x  10“^/?}  and  is  stretched  in 
the  normal  direction  from  an  initial  value  o’.  Ay,  -  355  x  lO'^m  to  a  final  value  of 

Ay/  =  5.594  X  10”^n?.  Grid  densities  are  ba  ,cd  upon  a.ssuming  a  Blasius  flat-plate 
boundary-layer  thickness  at  the  sampling  location  (49): 


v/7?h: 


(6.30) 


where  .r  is  the  distance  to  the  sampling  location  and  Rcx  is  the  Reynolds  number 
based  upon  the  velocity  behind  the  moving  shock  wave.  The  initial  spacing  in  the 
normal  direction  is  chosen  as  Ay,  »  6/20  .  The  axial  and  final  normal  spacing  are 
ch.o.'ten  to  provide  a.  reasonable  aspect  ratio  throughout  l.hc  domain. 

Solutions  i^resented  can  again  be  taken  as  the  results  obtained  through  uti¬ 
lizing  either  .ATNSCi  or  .ATNSC2  since  the  steady  Ileynolds  nuinbei  is  suinciently 
large  as  to  provide  for  the  exact  same  .solution  from  cithei  algorithm.  Figure  6.32 
is  the  solution  obtained  when  the  c  values  were  held  iIk'  .^aine  as  foi  the  shock- 
bou'!  .iary -layer  interaction  test  case  of  the  previous  section.  \'tt!iiely.  r,  =  (3  =  0.0 
and  £2  =  c>i  =  0.025.  The  numerical  solution  is  compared  to  the  theory  of  .Mirels  [31] 
whicli  is  valid  for  weak  shocks.  Time  is  referenced  to  the  time  of  the  shock  wave 
passing  the  sampling  point.  The  peak  heat  flux  is  much  less  than  theory  or  exper¬ 
iment.  but  it  should  be  realize,  that  the  theoretical  value  goes  to  infinity  at  the 
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exact  moment  of  shock  passage.  The  numerical  solution  is  initially  less  than  theory, 
but  soon  turns  so  as  to  become  greater  that  the  theoretical  flux.  The  numerical 
solution  continues  to  be  greater  than  the  theoretical  value  as  the  shock  continues  to 
progress  downstream.  Agreement  between  the  numerical  solution  and  experiment 
is  good  beyond  approximately  0.3  msec.  Figure  6.33  gives  the  percent  difference, 
Eq  3.17,  between  the  theoretical  value  and  the  .ATNSC  solution  as  a  function  of 
time.  The  initial  large  difference  is  expected  as  theory  predicts  an  infinite  heat  flux. 
The  numerical  solution  continues  to  exhibit  a  larger  percent  difference  as  time  goes 
on.  A  dramatic  change  in  the  slope  of  the  heat  flux  profile  at  approximately  O.Oo 
msec  suggests  there  may  be  an  anology  to  the  behavior  of  the  skin  friction  profile  of 
the  previous  sectioii  with  regards  to  c. 

Turning  again  to  the  work  of  Seider  and  Hanel  [39],  they  observed  the  wall 
temperature  on  a  flat  plate  approached  the  freestream  temperature,  rather  than 
the  recovery  temperature,  as  c  was  increased.  They  also  found  that  this  effect  was 
negated  when  they  increased  the  number  of  points  in  the  boundary-layer  from  7  to 
14.  While  increasing  the  number  of  points  in  the  boundary-layer  is  rea.sonable  foi 
steady  slate  calculations,  for  the  unsteady  calculations  of  this  test  case  there  is  no 
boundary-layer  immediately  at  shock  passing. 

Figures  6.34  and  6.35  depict  the  .ATW’SC  solution  when  c  is  set  to  0  for  all  ficULs. 
The  peak  heat  flux  has  increased  from  the  value  of  0.94  BTU/./7^  ■  s  of  Figure  6.32 
to  1.24  BTU//^^  •  The  heat  flux  profile  is  much  smoother  than  Figure  6.3=1  and 
approaches  the  theoretical  values  vcia  rapidlx.  FiqureC.35  shows  that  the  caknhiled 
heat  flux  becomes  greater  than  the  theoielic.il  value  atid  then  drops  to  a  value  ~'/i 
lower  that  theoretical,  followed  by  a  ri.se  to  a  final  level  just  1%  lower  that  the 
theoretical  level.  The  early  rise  above  theoretical  levels  followed  by  a  drop  below 
the  theoretical  level  may  perhaps  be  explained  by  examining  Figure  6.34.  The 
theoretical  curve  must  approach  infinity  ,i.s  time  goes  to  zero.  Following  backwards 
along  the  curve,  the  rise  in  heal  Iransfei  must  begin  before  the  associated  rise  in 
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Figure  6.32.  Meat  Flux  History  (ej  =  £3  =  O.ca  =  q  =  0.025) 


Figure  6.33.  Percent  Difference  Between  Theory  and  .ATN.SC  .Solution 
(c,  =  C3  =  0.  C2  =  c.|  =  0.025) 


103 


the  numerical  solution,  or  any  physical  solution,  since  neither  approach  infinilv  at 
time  zero.  Both  the  numerical  and  theoretical  heat  flu.\  profiles  are  shifted  below 
the  e.Kperimental  profile.  Increasing  the  shock  Mach  nuinbei  to  the  upper  end  of  the 
e.xperimental  tolerance  is  suggested,  and  will  be  examined  shortly. 

Visbal’s  implementation  of  the  Beam-Warming  algorithm,  Ecj  -3.39,  is  also  used 
to  obtain  a. solution  to  this  test  case.  Trapezoidal  time  differencing  is  used  to  provide 
second-order  accuracy  in  time. 

Figures  6.36  and  6.37  repre.sent  the  solution  obtained  using  the  iiominal  rec¬ 
ommended  value  of  the  second  and  fourtli  older  damping  coefficients  as  given  in 
.Section  6.3.  Boam-W^nrming  piedicU  a  peak  heat  flux  between  that  of  Figures  6.3! 
and  6.32.  Rather  severe  oscillations  occur  in  the  heat  flux  just  after  the  shock  passes 
the  sampling  point  but  eventually  damp  out  to  a  final  heat  flux  value  26%  higher  than 
the  theorectical  value.  Beam- Wanning  heat  flux  values  agree  well  with  experiment 
after  approximately  0.2  msec,  'his  is  surprising  in  light  of  the  severe  oscillations 
occuring  prior  to  this  time. 

The  values  of  tlie  damping  coefficient.*?  are  now  doubled,  consistent  with  the 
shock- boundaiA  layer  interaction  case,  and  the  solution  of  Figures  6.3S  and  6.39 
obtained.  'I’he  peak  heal  flux  has  been  reduced  .significantly  and  the  o.sfiUation.s 
damped  somewhat,  but  the  final  heat  flux  value  has  increased  to  33%-  above  the 
theoretical  value.  Beam  Warming  predictions  are  now  higher  than  the  experimc*ntai 
values  after  approximately  0.3  ms^c.  Thus,  as  the  oscillations  are  clamped  «iUt.  the 
overall  heat  flux  prediction  .suffers. 

The  clamping  coefficients  are  linalh  reduced  to  half  of  their  nominal  value 
resulting  in  the  solutions  of  Figures  6.10  and  6.11.  The  i)eak  heat  flux  is  apiiroaciiiiig 
a  value  consistent  with  the  results  of  Figure  6.34  but  the  oscillations  have  become 
particularly  .severe.  In  fact,  reducing  the  damping  much  below  thi.s  point  tau.'.es  ilu- 
algorithm  to  become  unstable.  'I'he  lower  damping  of  this  solution  provide.s  the  In^st 
Beam  Warming  prediction  of  the  final  In-ai  flu.x  value,  but  thi.s  i.s  stll  1N%  above  the 


Figure  6.36;  Heat  Flux  Hi.story  (w^  =  0.02,  =  0.04,  tuf  =  w”  =  0,25) 


Figure  6.37.  Percent  Difference  Between  Theory  and  Beam- Wanning  Solution 
■  (w<  =  0.02,  cf  =  0.04,  o-f  =  =  0.25) 
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theoretical  value.  Oscillations  in  the  Beam-Warming  solutions  of  Figures  6.36,  6.38, 
and  6.40  suggest  a  common  frequency,  independent  of  the  damping  applied.  Further 
examination  of  the  Beam- Warming  solutions  confirm  this.  Beam- Warming  bieaks 
the  initial  shock  wave  into  a  series  of  compressions  and  expansion  shocks  propagating 
downstream.  As  these  nonphysical  waves  pass  the  sampling  point,  oscillations  in  the 
heat  flux  are  observed.  Varying  the  amount  of  damping  applied  affects  the  magnitude 
of  the  jumps  across  these  waves,  but  not  the  frequency  at  which  they  are  generated. 

A  new  grid  is  now  utilized  in  an  attempt  to  see  the  effect  on  heat  flux  of 
decreasing  Ax  in  the  vicinity  of  the  sampling  point.  This  new  grid  is  shown  in 
Figure  6.42.  The  grid  spacing  varies  in  the  axial  direction  from  an  initial  value  of 
Ax  =  5.066  X  10“^  at  the  upstream  boundary  to  Ax  =  1.252  x  lO"'*  at  the  sampling 
point  to  Ax  =  2.026  x  10""  at  the  downstream  boundary.  The  grid  is  stretched  away 
from  the  wall  with  the  same  initial  spacing  as  the  grid  of  Figure  6.31  but  with  a  final 
spacing  of  Ay  =  1.252  x  lO"'^  at  the  upper  boundary. 

Figures  6.43  and  6.44  depict  I  he  result  of  utilizing  ATNSC  to  arrive  at  a  solution 
with  e  =  0  in  all  fields.  Figure  6.43  shows  that  the  predicted  peak  heat  flux  has 
increased  24%  over  that  of  Figure  6.34.  The  heat  flux  profile  seems  to  have  been 
stretched  immediately  at  shock  passage,  lemaining  unchanged  after  appioximatelx 
0.1  msec,  and  is  still  shifted  below  the  experimental  profile.  This  is  verified  by 
comparison  of  Figure  6.44  and  Figure  6.35. 

The  Beam- Warming  algorithm  is  applied  to  this  new  grid  using  the  nominal 
damping  coefficients  mentioned  earlier.  Results  are  shown  in  Figures  6.-15  and  6.16. 
While  th'’  peak  predicted  heat  flux  has  increased  over  that  of  Figure  6.36,  it  is  just 
beginning  to  approached  the  level  of  the  .-\TNSC  prediction  on  the  previous  grid.  At 
times  greater  than  0.3  msec,  the  solution  is  tending  awav  form  the  good  agreement 
with  experimental  data  exhibited  in  Figure  6.36.  Figure  6.46  shows  a  8%  inciea.sc  in 
the  difference  from  theory,  as  compared  to  Figure  6.37,  for  the  final  heat  flux  \alne. 
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Figure  6.42.  Second  Grid  for  Heat  Flux  Solutions 


Finally,  solutions  were  computed  with  the  .ATNSC  and  Beam-Warming  algo¬ 
rithms  with  Ms  =  1.117  for  comparison  against  each  other  and  the  experimental 
data.  ATNSC  calculations  again  used  e  =  0,  while  nominal  values  of  the  damping 
coefficients  weie  used  for  the  Beam-Waiming  solution.  Figure  6.47  presents  the  time 
history  for  this  case.  Agreement  between  theory,  experiment,  and  the  .ATNSC  solu¬ 
tion,  is  excellent  in  this  case,  and  suggests  that  the  experimental  shock  Mach  numbei 
was  probably  closer  to  2%  above  the  nominal  value.  Beam-Warming  again  undei pre¬ 
dicts  the  peak  heat  flux,  displays  oscillations  as  compre.ssions  and  expan.sion  shock.s 
pass  the  sampling  location,  and  utimately  overpredicts  the  heat  flux  level. 

It  is  clear  that  a  TVD  based  algorithm  such  as  .ATNSC  is  necessary  to  obtain 
acceptable  solutions  for  the  problem  of  unsteady  shock-induced  heat  transfer.  The 
solutions  are  free  of  any  oscillations  and  come  extremely  close  to  the  theoretical 
values  of  Mirels  [30].  In  contrast,  the  Beam- Warming  algoi  ithm  yields  solutions  with 
oscillatory  behavior,  due  to  the  generation  of  nonphysical  waves.  e\en  at  ielati\el\ 
high  values  of  the  damping  coefficients.  In  no  instance  did  the  Beam- Warming 
algorithm  yield  an  acceptable  comparison  with  theoretical  values,  predicting  heat 
.Hux  values  IS  —  34%  higher  than  theory. 


.All  AT.NSC  solutions  are  undertaken  at  a  CFL  of  0.96.  with  the  time  .step 
restriction  of  Eci  1.2S.  This  results  in  appro.viaiately  7000  sweeps  to  arrive  a(  a  lime 


Figure  6.45.  Heat  Flux  History  (wf  =  0.02,  w 
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Figure  6.47.  Heat  Flux  History  {Mg  =  1.117) 


of  1.0  msec  for  the  first  grid,  and  approximately  14000  sweeps  for  the  second  grid. 
The  Beam- Warming  solutions  utilize  a  constant  time  step  which  was  chosen  from 
the  smallest  time  step  taken  by  the  .ATN.SC  algorithm. 

6.0  Conclusions 

TVD  methodolgy  has  been  applied  to  problems  not  previously  examined  using 
TVD  schemes.  Prior  research  concenlraled  mainly  on  supersonic  and  hypersonic 
flows,  both  inviscid  and  viscous,  and  was  almost  solely  directed  towaid  obtaining 
steady-state  solutions.  The  current  effort  has  extended  the  T\‘D  mcthodologx  to 
inviscid  transonic  ca.scade  flows,  \iscous  flows  with  shock  induced  laminar  boundaiy 
layer  separation,  and  unsteady  laniinai  flows  wdth  signifeant  shock-induced  heat 
transfer.  .Additionally,  an  algorithm  was  developed  herein  that  shows  promise  for 
application  to  low  Reynolds  number  situations. 

Transonic  cascade  flows  arc  currently  of  great  interest  to  gas  turbine  engine 
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designers  and  researchers  [2,  14.  13).  Analysis  of  these  flows  is  a  severe  test  of  an 
algorithm  because  of  the  wide  .Mach  number  range,  typically  0.3  <  M  <  1.3  ,  and 
the  fact  that  the  flow  is  confined  in  a  passage  where  wave  systems  tend  to  reflect  back 
into  the  domain.  Results  presented  in  Section  3.4  show  that  not  only  do  the  TVD 
schemes  of  Sections  2.2.2  and  2.2.3  yield  steady-state  results  in  e.xcelleni  agreement 
with  experiment,  but  that  the  transient  behavior  is  also  modeled  correctly.  Correct 
modeling  of  this  transient  behavior  is  an  important  achievement,  since  previous 
efforts  have  been  unsuccessful  in  this  regard  (38,  15). 

Laminar  shock-boundary-layer  interaction  has  been  studied  by  numerous  re¬ 
searchers  using  highly  regarded  algorithms  such  a.s  the  .MacCormack  [28],  Dawes  [10], 
Beam- Warming  [11],  and  Newton  [2G|  methods.  While  acceptable  predictions  of 
the  pressure  profile  in  the  boundary-layer  have  been  coi..puted,  researchers  have 
remained  unable  to  accurately  compute  the  skin  friction  profile.  The  ATNSC  algo¬ 
rithms  finally  provide  the  means  to  accurately  compute  pressure  profiles,  separation 
and  reattachment  locations,  and  skin  friction  piofiles  in  excellent  agreement  with  the 
available  experimental  data.  Results  given  in  Section  6.3  are  testament  to  this.  In 
addition  to  the  success  in  solving  this  complex  problem,  the  current  effort  extends 
the  knowledge  of  how  TVD  entropy  coijcclion  affects  the  boundary-layer  velocity 
profile. 

Unsteady  shock-induced  heat  transfer  has  been  studied  theoretically  (30,  31), 
experimentally  [1-1.  41),  and  has  rccenth  become  of  interest  computationally.  The 
increa.sed  interest  is  due  to  the  enhanced  system  performance  acailable  through  ac 
curate  knowlege  of  the  heal  transfer  (G).  However,  it  is  not  uncommon  for  com 
puled  heat  flux  values  to  be  an  order  of  magnitude  different  from  experimental  val¬ 
ues  (17).  The  .M’NSC  algorithms  represent  a  significant  advancement  in  the  state- 
of-the-art  for  computing  shock  induced  Inut  transfer.  .Solutions  computed  with  the 
.ATN.SC  schemes,  presented  in  .Section  6.4.2.  are  far  superior  to  the  nonphysical 
Beam -Warming  solulions  and  agree  well  with  both  theory  and  experiment.  .Addi 


Uonally.  the  ATNSCl  data  processing  rate  is  16%  faster  than  that  for  Bearn-Warming 
using  the  same  time  step. 

Results  of  applying  the  ATN!:)C2  scheme  to  the  viscous  Burgers'  equation  sug¬ 
gest  that  the  scheme  may  be  beneficial  to  investigators  studying  low  Reynolds  num¬ 
ber  flows.  The  scheme  is  extremely  stable  and  exhibits  superb  accuracy  for  this  model 
problem.  Conventional  wisdom  states  that  TVD  schemes  are  not  recommended  for 
problems  containing  no  discontinuities,  since  they  degrade  to  first-  order  accuracy- 
near  points  of  e.xtrema  [4.5].  Behavior  of  the  ATNSC2  scheme  suggests  that  TVD 
algorithms  warrant  a  closer  examination  for  application  to  shock-free  flows. 

The  research  described  herein  represents  a  significant  contribution  to  the  field 
of  computational  fluid  dynamics.  Suggestions  for  extending  the  research  efforts  are 
presented  in  the  following  section. 

6.6  Furllier  Research 

Further  research  needs  to  be  conducted  in  several  areas.  The  main  (|ue.‘«tion 
in  terms  of  the  inviscid  investigations  relates  to  the  effect  of  grid  cell  .skcwne.ss  on 
entropy  production.  Isolating  this  investigation  to  the  inviscid  ca.se  will  enable  any 
insights  to  be  directly  implemented  into  •.  i  Iscous  algorithm,  were  it  would  be  difficult 
10  distinguish  between  spurious  entropy  production  and  normal  viscous  dissipation. 

Two  areas  should  be  addre.ssed  for  the  viscous  algorithms.  First  is  th'^  effect 
of  the  second-order  accurate  algorithm  at  low  Reynolds  numbers.  'Phe  Reynolds 
munijers  of  interest  herein,  except  for  Burgers'  et|uation,  were  on  the  ordei  of  lO"* 
and  no  distinguishable  difference  between  the  ATN.SC1  and  .AT.\'.SC2  algorithms  was 
observed,  llow'ever,  the  low  Reynolds  number  test  provided  by  Burgeis'  eciuation 
dramatically  showed  the  enhanced  performance  of  the  second -order  algorithm  over  its 
first  order  counterpart.  .Although  no  low  Reynolds  number  flows  conlaiuiug  shock.s 
comes  to  mind,  except  for  the  high  altitude  transistion  regime  were  the  continuum 
assumption  fails  to  hold,  further  investigation  into  the  applicability  of  the  .srcoml- 


order  algorithm  is  certainly  warranted.  Further  analysis  of  the  truncation  error 
cancelling  terms  in  ATNSC2  should  be  performed.  If  some  of  the  terms  are  negligible 
for  certain  flow  situations,  a  dramatic  increase  in  computational  efficiency  maj  be 
possible. 

Second,  performance  of  the  .ATNSC  algorithms  with  a  turbulence  model  in¬ 
cluded  should  be  investigated.  .At  higher  overall  pressure  ratios  than  the  one  con¬ 
sidered  in  Section  6.3.  the  boundary-layer  tends  to  transition  to  turbulent  upon 
reattachment,  resulting  in  higher  skin  friction  levels.  As  previously  mentioned,  non- 
TVD  algorithms  under-  predict  the  laminar  skin  friction  level  beyond  reattachnicnt. 
Preliminary  investigations  show  that  implementation  of  a  Baldwin-Loma.x  turbu¬ 
lence  model  [3]  tends  to  ovcrpredict  the  skin  friction  in  the  turbulent  region.  This 
initially  leads  the  author  to  believe  that  turbulence  models,  developed  in  the  con- 
te.\t  of  non-TVD  algorithms,  may  compensate  for  a  non-TVD  scheme’s  tendency  to 
under-predict  skin  friction. 

Overall,  the  TVD  based  viicous  algorithms  perform  e.Kceptionally  well  on  the 
test  cases  herein.  Emphasis  must  be  placed  on  applying  the.se  algorithms  to  even 
more  rigorous  lest  cases  so  as  to  gain  an  even  greater  understanding  of  theii  weak¬ 
nesses  as  well  as  their  strengths. 
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Appendix  A.  VtscoiLs  Jacobians 


A, (LI)  A,(i,2)  /L(L3)  /L(L.'l) 
A, (2,1)  -4.(2, 2)  A, (2, 3)  A, (2,^1) 
Ax(3,i)  A, (-3, 2)  A, (3, -3)  Ai(3,4) 
Ai(.Ll)  A, (4, 2)  A, (4, 3)  A, (4, 4) 


>1.(1,  i)=o 


A.(2, 1)  =  ^  {f)  ^(2/.  a-  A)  +  I;  (a)  g  -  (2/;  -1-  A)^  (^)  -  . 

>1.(3, 1}  =  S  (7)  +  fe  (^)l  -  [i  if)  +  I  if)] 

A, (4, 1)  =  f  A.(2, 1)  -  f  [(2/.  +  X)£  if)  +  A|;  (^)] 

(f)  -  5^  (^)l  -  Mi  (^)  -  i  (^)l} 
(I  -  M  (i  (7)  i  {7)1  -  >  (i  (?)  i  (?^)1 

A.(l,2)=0 

•■'.(2,2)  =  ^  (f)  |:(2,.  A)  -f  i  (";)  |4  +  (2,,  ,\)l  (1) 

•■'.(■5,2)  =  fete  (:)  +  !;  (7)1 +<‘l;(i) 

.4, (4, 2)  =  54, (2, 2)  .f  i  [(2,.  A);|  (f)  .f  Ai  (=)] 

+i/c4fe  te  (5)  -  ji  {“^)1  - (?)} 


1,(1,3)=0 


'.(•2.3)  =  fete  (;)  +  !;  (7)1 -K-i 


A>>  ™ 


OF,. 

dUr 


(A.2) 


is  given  by 


-(2/^  +  A)^ 


P 


0 

0 


0 

0 

0 


-'13  = 


dU„ 


-A^ 


-(/i  +  A)^ 


0  0  0 


0 

p 


0 

0 


„iL  \”l  n 

/'  fi!  •' 


(A.:J) 


A-II.S 


Bi  = 


dG, 

dU 


/?.(!,  1)  /i.(L2)  i?,(L3)  5, (1,4) 
A(2,l)  5,(2, 2)  5, (2, 3)  5, (2, 4) 
5,(3, 1)  5,(3,2)  5,(3,3)  5, (3,4) 
5,(4, 1)  5, (4, 2)  5,(4, 3)  5, (4,4) 


(A.4) 


5,(1, 1)  =  0 


E  [a  fa)  4.  |L 

p  11'-  \p/  t'y 

(7)1-"!A(?)+A(?)1 

L  (la^  ^  '  J. 

X  \  P  )  9p  Oy 

(?)A(2/-  +  -')-4(p)-(2"  +  '')| 

0.(3,!)- A1 

(2"+4(?)+4(7)i 

{S  [A  (f )  -  5  A  (^)l  -  [^,  if)  -  ^  (^)l } 

^  ^  (  hM  _  Eltf  f-2.  ^  JS.)  •  ®  Vl 

•  p  p)  [9x\p)  ■  Vp;j  l9rU*/  •  9y\p-)\ 

Bi(i,2)  =  0 

».(2.2)  =  £  [A  (;)  +  !(?))+ 4  (i) 

»■  2)  =  A  (?)  A(2'‘  --')  +  *  (?)  t  +  4  (?) 

0,01/2)  =  ;0.(3,2)  +  (J  +  [A  (;)  +  A  (^)l 

-i-'/f.-  {s  (A  (?)  -  ? A  (^)l  -  4  (?)}  +  "75  (?) 

5,(1,3)  =  0 

e.(2-3)  =  i(A(?)  +  A(7)j  +  "A(?) 

0.(3,3)  =  A  (?)  A(2"  +  -M  +  5  (?)  5  +  (2"  ■!■  -'jA  (?) 

0.01,3)  =  JB,i3,3)  +  1  [(2„  +  A)A  {;)  +  aA  (?)! 

+‘/'^» {i  lA  (?)  -  4  (=^)l  - *-A  (?)} 

+7{i?lA(?)  +  A(7)l+4(?)} 

5,(1,4)  =  0 


5,  (2. 4) 


A- 11 5) 


5i(3,4)  = 
B,(4,4)  = 


is  given  by 


-  (2/(  +  A  -  ^ 


B2 


dG, 

dUy 


(A.5) 


0 


-(2/4 +A)^ 


Jls. 

C,  pi 


0  0  0 

0  0 
p 

0  0 
~  ^  +  A  -  ^  ^  _ 


0  0  0  0 

0  f  0 

-A^  A  0  0 

-in  +  xy-f  Xf  fx^  0 


(A.6) 
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Appendix  B.  ATEC  Routines  and  FLOWTRACE  Results 


B.l  Description  of  ATEC  Routines 

For  clarity,  the  routines  are  described  in  the  order  they  would  generally  be 
called,  independent  of  the  sweep  direction. 

ATEC  -  main  program 

GEOMETRY  -  computes  cell  centers  based  on  corner  values 
TFORM  -  computes  the  metric  transformation  terms 
INITIAL  -  enforces  the  initial  conditions 

STORE  -  stores  the  dependent  variables  at  the  current  time  level 
FLUXF  -  computes  the  direction  flux 

FLUXG  -  computes  the  ?/  direction  flux 

ROEAVGZ  -  computes  Roe  averaged  quantities  along  constant  //  line.s 
ROEAVGE  -  computes  Roe  averaged  quantities  along  constant  ^  lines 
EVALUEZ  -  computes  the  ^  eigenvalues 
EVALUEE  -  computes  the  ■/;  eigenvalues 
TMSTEP  -  computes  the  allowable  time  step 

ALPHAZ  -  computes  the  difference  of  characteristic  variables  in  the  f  direction 

ALPH.AE  -  computes  the  difference  of  characteristic  variables  in  the  ;/  direction 

GCALCZ  -  computes  the  flux  limiters  for  the  f  direction  flux 

GC.ALCE  -  computes  the  flux  limiters  for  the  t)  direction  flux 

BET.ATVDZ  -  computes  artificial  dissipation  for  ^  sweep 

BETATVDE  -  computes  artificial  dissipation  for  sweep 

EVECTORZ  -  computes  the  eigenvectors  for  the  E,  eigenvalues 

EVECTORE  -  computes  the  eigenvectors  for  the  p  eigenvalues 

ARTCOMPZ  -  computes  the  final  artificial  dissipation  for  the  £,  direction  sweep 

.ARTCOMPE  -  computes  the  final  artificial  dissipation  for  the  //  direction  sweep 


FSOLVE 

GSOLVE 

BNDBLD 

BNDEX 

BNDPER 

BNDIN 

NORM 

OUTPUT 


-  solves  for  the  drpenclent  variables  during  the  ^  sweep 

-  solves  for  the  dependent  variables  during  the  //  sweep 

-  enforces  the  blade  or  wall  surface  boundary  conditions 

-  enforces  the  exit  plane  boundary  conditions 

-  enforces  the  periodic  boundary  conditions 

-  enforces  the  inlet  plane  boundary  conditions 

-  computes  the  L-2  and  norms 

-  outputs  the  solution  vector 
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B.2  ATEC-FV  ( Ft nUe- Volume  Formulation) 

The  data  processing  rate  tor  the  finite- volume  formulation  is  1.2274  x  10"® 
seconds  per  grid  point  per  time  level  for  the  CRAY  X-MP/216,  utiiling  a  177  x  20 
grid.  FLOWTRACE  results  are  for  1000  iterations  (2000  time  levels). 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  TIME  USED  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

Ac  cum'/, 

BETATVDE 

1.40E+01 

352000 

3.98E-05 

16.10 

16.10 

BETATVDZ 

1.15E+01 

57000 

2.02E-04 

13.23 

29.33 

GSOLVE 

6.99E+00 

352000 

1.99E-05 

8.04 

37.37 

ATEC 

6.55E+00 

1 

6.55E+00 

7.54 

44.91 

ROEAVGE 

5.55E+00 

3000 

1.85E-03 

6.39 

51.30 

ROEAVGZ 

5.30E+00 

3000 

1.77E-03 

6.10 

57.40 

FSOLVE 

4.86E+00 

57000 

8.52E-05 

5.59 

62.99 

EVECTORE 

3.52E+00 

352000 

l.OOE-05 

4.05 

67.04 

ALPHAZ 

3.28E+00 

3000 

l.lOE-03 

3.78 

70.82 

GCALCZ 

3.16E+00 

3000 

1.05E-03 

3.63 

74.45 

ARTCOMPE 

3.01E+00 

352000 

8.55E-06 

3.46 

77.91 

EVECTORZ 

2.38E+CO 

57000 

4.18E-05 

2.74 

80.65 

FLUXG 

2.37E+00 

5000 

4.74E-04 

2.72 

83.37 

ARTCOMPZ 

2.32E+00 

57000 

4.06E-05 

2.66 

86.04 

ALPHAE 

2.29E+00 

2000 

1 . 14E-03 

2.63 

88.67 

FLUXF 

2.27E+00 

5000 

4.54E-04 

2.61 

91.28 

GCALCE 

2.15E+00 

2000 

1.07E-03 

2.47 

93.75 

NORM 

1.06E+00 

400 

2.65E-03 

1.22 

94.97 
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TMSTEP 

.  1. 

03E+00 

IGOO 

1. 

,03E- 

•03 

1. 

19 

96. 

15 

OUTPUT 

7, 

.18E-01 

1 

7. 

.18E- 

■01 

0, 

,83 

96. 

,98 

EVALUEE 

6. 

,50E-01 

3000 

2, 

.  17E- 

•04 

0. 

,75 

97. 

,73 

evaLuez 

6. 

,18E-01 

3000 

2, 

.06E- 

•04 

0, 

,71 

98, 

.44 

BNDBLD 

4. 

,66E-01 

5000 

9. 

.32E- 

■05 

0, 

,54 

98. 

.97 

BNDEX 

3. 

,91E-01 

5000 

7, 

.82E- 

-05 

0, 

,45 

99. 

.42 

BNDPER 

3, 

.35E-01 

5000 

6, 

.70E- 

■05 

0, 

.39 

99, 

.81 

BNDIN 

1, 

.46E-01 

5000 

2 

.92E- 

-05 

0, 

.17 

99, 

.98 

STORE 

1, 

.61E-02 

100 

1 

.61E- 

-04 

0 

.02 

100, 

.00 

TFORM 

3, 

.25E-03 

1 

3 

.25E- 

-03 

0 

.00 

100, 

.00 

GEOMETRY 

2 

.75E-04 

1 

2 

.75E- 

-04 

0 

.00 

100 

.00 

INITIAL 

1 

.83E-04 

1 

1 

.83E- 

-04 

0 

.00 

100 

.00 

Totals 

8 

.69E+01 

1689505 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  'IN-LINE’  FACTOR  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 

(Factors  Greater  Than  1  Could  Indicate  Candidates  for  In-Lining) 


Routine  Name 

Tot  Time 

#  Calls 

Avg  Time  Percentage 

"In-Line" 

ARTCOMPE 

3.01E-<-00 

352000 

8.55E-06 

3.46 

349.66 

EVECTORE 

3.52E-<-00 

352000 

l.OOE-05 

4.05 

298.90 

GSOLVE 

6.99E-t-00 

352000 

1.99E-05 

8.04 

150.58 

BETATVDE 

1.40E■^01 

352000 

3.98E-05 

16.10 

75.20 

ARTCOMPZ 

2.32E-<-00 

57000 

4.06E-05 

2.66 

11.92 

EVECTORZ 

2.38E-»-00 

57000 

4.18E-05 

2.74 

11.59 

B-124 


FSOLVE 

BETATVDZ 

BNDIN 

4.86E+00 

1.15E+01 

1.46E-0i 

57000 

57000 

5000 

8.52E-05 

2.02E-04 

2.92E-05 

5.59 

13.23 

0.17 

5.68 

2.40 

1.45 

BNDPER 

3.35E-01 

5000 

6.70E-05 

0.39 

0.63 

BNDEX 

3.91E-01 

5000 

7.82E-05 

0.45 

0.54 

BNDBLD 

4.66E-01 

5000 

9.32E-05 

0.54 

0.46 

EVALUEZ 

6.18E-01 

3000 

2.06E-04 

0.71 

0.12 

EVALUEE 

6.50E-01 

3000 

2.17E-04 

0.75 

0.12 

FLUXF 

2.27E+00 

5000 

4.54E-04 

2.61 

0.09 

FLUXG 

2,37E+00 

5000 

4.74E-04 

2.72 

0.09 

GCALCZ 

3.16E+00 

3000 

1.05E-03 

3.63 

0.02 

ALPHAZ 

3.28E+00 

3000 

l.lOE-03 

3.78 

0.02 

GCALCE 

2.15E+00 

2000 

1.07E-03 

2.47 

0.02 

ALPHAE 

2.29E+00 

2000 

1 . 14E-03 

2.63 

0.01 

ROEAVGZ 

5.30E+00 

3000 

1.77E-03 

6.10 

0.01 

ROEAVGE 

5.55E+00 

3000 

1.85E-03 

6.39 

0.01 

ATEC 

6.55E+00 

1 

6.55E+00 

7.54 

0.00 

NORM 

1.06E+00 

400 

2.65E-03 

1.22 

0.00 

TMSTEP 

1.03E+00 

1000 

1.03E-03 

1.19 

0.01 

OUTPUT 

7.18E-01 

1 

7.18E-01 

0.83 

0.00 

STORE 

1.61E-02 

100 

1.61E-04 

0.02 

0.01 

TFORM 

3.25E-03 

1 

3.25E-03 

0.00 

0.00 

GEOMETRY 

2.75E-04 

1 

2.75E-04 

0.00 

0.00 

INITIAL 

1.83E-04 

1 

1.83E-04 

0.00 

0.00 

Totals 

8.69E+01 

1689505 
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B.3  ATEC-FD  (Chain-Rule  Formulation) 

The  data  processing  rate  for  the  chain-rule  formulation  is  1.24 15  x  10“’  seconds 
per  grid  point  per  time  level  for  the  CRAY  X-MP/216,  utiiling  a  177  x  20  grid. 
FLOWTRACE  results  are  for  1000  iterations  (2000  time  levels). 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  TIME  USED  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 


Routine  Name 

Tot  Time 

#  Calls 

Avg  Time  Percentage 

Accum'/, 

BETATVDE 

1.37E+01 

352000 

3.88E-05 

15.56 

15.56 

BETATVDZ 

1.13E+01 

57000 

1.97E-04 

12.81 

28.37 

GSOLVE 

8.01E+00 

352000 

2.28E-05 

9.12 

37.49 

ATEC 

6.58E+00 

1 

6.58E+00 

7.49 

44.99 

FSOLVE 

5.59E+00 

57000 

9.80E-05 

6.36 

51.34 

ROEAVGE 

5.54E+00 

3000 

1.85E-03 

6.30 

57.64 

ROEAVGZ 

5.29E+00 

3000 

1.76E-03 

6.02 

63.67 

EVECTORE 

3.54E+00 

352000 

l.OlE-05 

4.03 

67.70 

ALPHAZ 

3.25E+00 

3000 

1.08E-03 

3.69 

71 .39 

GCALCZ 

3.15E+00 

3000 

1.05E-03 

3.58 

74.97 

ARTCOMPE 

2.99E+00 

352000 

8.49E-06 

3.40 

78.37 

FLUXG 

2.36E+00 

5000 

4.72E-04 

2.69 

81.06 

EVECTORZ 

2.34E+00 

57000 

4.11E-05 

2.67 

83.73 

ARTCOMPZ 

2.26E+00 

57000 

3.97E-05 

2.57 

86.30 

FLUXF 

2.26E+00 

5000 

4.52E-04 

2.57 

88.88 

ALPHAE 

2.26E+00 

2000 

1.13E-03 

2.57 

91.45 

GCALCE 

2.14E+00 

2000 

1.07E-03 

2.44 

93.89 

NORM 

1.03E+00 

400 

2.57E-03 

1.17 

95.06 
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TMSTEP 

1. 

.03E+00 

1000 

1, 

.03E-03 

1, 

.17 

96. 

.23 

OUTPUT 

7. 

,13E-01 

1 

7. 

,13E-01 

0, 

.81 

97, 

.04 

EVALUEE 

6. 

.33E-01 

3000 

2, 

.llE-04 

0, 

.72 

97. 

.76 

EVALUEZ 

6. 

.15E-01 

3000 

2. 

.05E-04 

0, 

.70 

98. 

.46 

BNDBLD 

4. 

.63E-01 

5000 

9, 

,26E-05 

0, 

.53 

98. 

.99 

BNDEX 

3, 

.91E-01 

5000 

7, 

,82E-05 

0, 

.44 

99, 

.43 

BNDPER 

3. 

.34E-01 

5000 

6, 

.68E-05 

0. 

.38 

99. 

.81 

BNDIN 

1, 

.46E-01 

5000 

2, 

.92E-05 

0, 

.17 

99, 

.98 

STORE 

1, 

.56E-02 

100 

1, 

.56E-04 

0, 

.02 

100, 

.00 

TFORM 

1 

.87E-03 

1 

1 

.87E-03 

0, 

.00 

100 

.00 

GEOMETRY 

2 

.79E-04 

1 

2 

.79E-04 

0, 

.00 

100, 

.00 

INITIAL 

2 

.05E-04 

1 

2 

.05E-04 

0 

.00 

100 

.00 

Totals  8.79E+01  1689505 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  'IN-LINE'  FACTOR  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 

(Factors  Greater  Than  1  Could  Indicate  Candidates  for  In-Lining) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

"In-Line" 

ARTCOMPE 

2.99E+00 

352000 

8.49E-06 

3.40 

352.44 

EVECTORE 

3.54E+00 

352000 

l.OlE-05 

4.03 

297.04 

GSOLVE 

8.01E+00 

352000 

2.28E-05 

9.12 

131.35 

BETATVDE 

1.37E+01 

352000 

3.88E-05 

15.56 

76.98 

ARTCOMPZ 

2.26E+00 

57000 

3.97E-05 

2.57 

12.20 

EVECTORZ 

2.34E+00 

57000 

4.11E-05 

2.67 

11.77 
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FSOLVE 

5.59E+00 

57000  9.80E-05 

6.36 

4.94 

BETATVDZ 

1.13E+01 

57000  1.97E-04 

12.81 

2.45 

BNDIN 

1.46E-0i 

5000  2.92E-05 

0.17 

1.46 

BNDPER 

3.34E-01 

5000  6.68E-05 

0.38 

0.64 

BNDEX 

3.91E-01 

5000  7.82E-05 

0.44 

0.54 

BNDBLD 

4.63E-01 

5000  9.26E-05 

0.53 

0.46 

EVALUEZ 

6.15E-01 

3000  2.05E-04 

0.70 

0.12 

EVALUEE 

6.33E-01 

3000  2.11E-04 

0.72 

0.12 

FLUXF 

2.26E+00 

5000  4.52E-04 

2.57 

0.09 

FLUXG 

2.36E+00 

5000  4.72E-04 

2.69 

0.09 

GCALCZ 

3.15E+00 

3000  1.05E-03 

3.58 

0.02 

ALPHAZ 

3.25E+00 

3000  1.08E-03 

3.69 

0.02 

GCALCE 

2 . 14E+00 

2000  1.07E-03 

2.44 

0.02 

ALPHAE 

2.26E+00 

2000  1.13E-03 

2.57 

0.02 

ROEAVGZ 

5.29E+00 

3000  1.76E-03 

6.02 

0.01 

ROEAVGE 

5.54E+00 

3000  1.85E-03 

6.30 

0.01 

ATECFD 

6.58E+00 

1  6.58E+00 

7.49 

0.00 

NORM 

1.03E+00 

400  2.57E-03 

1.17 

0.00 

TMSTEP 

1.03E+00 

1000  1.03E-03 

1.17 

0.01 

OUTPUT 

7.13E-01 

1  7.13E-01 

0.81 

0.00 

STORE 

1.56E-02 

100  1.56E-04 

0.02 

0.01 

TFORM 

1.87E-03 

1  1.87E-03 

0.00 

0.00 

GEOMETRY 

2.79E-04 

1  2.79E-04 

0.00 

0.00 

INITIAL 

2.05E-04 

1  2.05E-04 

0.00 

0.00 

Totals 

8.79E+01 

1689505 
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Appendix  C.  ATNSC  Routines  and  FLOWT'RACE  Results 

C.l  Description  of  ATNSCl  and  ATNSC2  Routines 

For  clarity,  the  routines  are  described  in  the  order  they  would  generally  be 
called,  independent  of  the  sweep  direction. 


.ATNSCl  (ATNSC2) 

GEOMETRY 

TEORM 

INITIAL- 

STORE 

EULFLUX 

YISFLUX 

ROE.AVGZ 

ROEAVGE 

EVALUEZ 

EVALUEE 

TMSTEP 

ALPHAZ 

ALPHAE 

GCALCZ 

GCALCE 

•JACOBIAN 

PSIZETA 

PSIETA 

BETATVDZ 

BETATVDE 

EVECTORZ 


-  main  program 

-  computes  cell  centers  based  on  corner  values 

-  computes  the  metric  transformation  terms 

-  enforces  the  initial  coiuliliuiis 

-  stores  the  dependent  variables  at  the  current  time  level 

-  computes  the  inviscid  flux  terms 

-  computes  the  viscous  flux  terms 

-  computes  Roe  averaged  quantities  along  constant  7/  lines 

-  computes  Roe  a.veraged  quantities  along  constant  ^  lines 

-  computes  the  ^  eigenvalues 

-  computes  the  ;/  eigenvalues 

-  computes  the  allowable  time  step 

-  computc.s  the  difference  of  characteristic  variables  in  the  i  direction 

-  compute.s  the  difference  of  characteristic  variables  in  the  //  direction 

-  computes  the  flux  limiters  for  the  ^  direction  flux 

-  computes  the  flux  limiters  for  the  ;/  direction  flux 

-  computes  the  vi.scous  -Jacobians  (ATNSC2  only) 

-  computes  the  final  viscous  flux  in  the  ^  direction 

-  computes  the  final  viscous  flux  in  the  r/  direction 

-  computes  artificial  dissipation  for  ^  sweep 

-  computes  artificial  dissipation  for  //  sweep 

-  computes  the  eigenvectors  for  the  ^  eigenvalues 
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EVECTORE 

ARTCOMPZ 

ARTCOMPE 

FSOLVE 

GSOLVE 

BNDBLD 

BNDEX 

BNDPER 

BNDIN 

NORM 

OUTPUT 


-  computes  the  eigenvectors  for  the  j/  eigenvalues 

-  computes  the  final  artificial  dissipation  for  the  ^  direction  sweep 

-  computes  the  final  artificial  dissipation  for  the  r/  direction  sweep 

-  solves  for  the  drpendent  variables  during  the  ^  sweep 

-  solves  for  the  dependent  variables  during  the  ;;  sweep 

-  enforces  the  blade  or  wall  surface  boundary  conditions 

-  enforces  the  exit  plane  boundary  conditions 

-  enforces  the  periodic  boundary  conditions 

-  enforces  the  inlet  plane  boundary  conditions 

-  computes  the  and  L,^  norms 

-  outputs  the  solution  vector 
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C.2  ATNSCl  (Constant  Damping) 

The  data  processing  rate  for  the  constant  c  case  is  1.6071  x  10“’  seconds 
per  grid  point  per  time  level  for  the  CR.AY  X-MP/216.  utiiling  a  133  x  60  grid. 
FLOWTRACE  results  are  for  200  iterations  (400  time  levels). 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  TIME  USED  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

Accum'/, 

VISFLUX 

7 .66E+00 

1000 

7.66E-03 

14.92 

14.92 

FSOLVE 

6.13E+00 

34800 

1.76E-04 

11.94 

26.86 

BETATVDZ 

4.69E+00 

34800 

1.35E-04 

9.13 

35.99 

GSOLVE 

4.40E+00 

52400 

8.40E-05 

8.57 

44.56 

BETATVDE 

4.25E+00 

52400 

8.11E-05 

8.28 

52.84 

ROEAVGE 

2.67E+00 

600 

4.45E-03 

5.20 

58.04 

ROEAVGZ 

2.66E+00 

600 

4.44E-03 

5.19 

63.23 

ATNSCl 

1.92E+00 

1 

1.92E+00 

3.74 

66.97 

ALPHAZ 

1.61E+00 

600 

2.68E-03 

3.13 

70.09 

EULFLUX 

1.59E+00 

1000 

1.59E-03 

3.11 

73.20 

GCALCZ 

1.53E+00 

600 

2.55E-03 

2.98 

76.18 

OUTPUT 

1.49E+00 

1 

1.49E+00 

2.90 

79.08 

PSIZETA 

1.20E+00 

600 

2.00E-03 

2.34 

81.42 

EVECTORZ 

1 . 19E+00 

34800 

3.41E-05 

2.31 

83.73 

ARTCOMPZ 

1 . 18E+00 

34800 

3.38E-05 

2.29 

86.02 

ALPHAE 

1.08E+00 

400 

2.71E-03 

2.11 

88.13 

THSTEP 

1.08E+00 

200 

5.38E-03 

2.09 

90.22 

GCALCE 

1.02E+00 

400 

2.55E-03 

1.99 

92.21 
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EVECTORE 

9.69E-01 

52400  1.85E-05 

1.89 

94.10 

ARTCOMPE 

8.09E-01 

52400  1.54E-05 

1.58 

95.67 

PSIETA 

8.05E-01 

400  2.01E-03 

1.57 

97.24 

NORM 

5.01E-01 

80  6.26E-03 

0.98 

98.22 

EVALUEZ 

3.24E-01 

600  5.39E-04 

0.63 

98.85 

EVALUEE 

3.23E-01 

600  5.39E-04 

0.63 

99.48 

BNDBLD 

1.97E-01 

1000  1.97E-04 

0.38 

99.86 

BNDEX 

4.17E-02 

1000  4.17E-05 

0.08 

99.94 

-BNDIN 

1.32E-02 

1000  1.32E-05 

0.03 

99.97 

STORE 

7.56E-03 

20  3.78E-04 

0.01 

99.98 

INITIAL 

5.34E-03 

1  5.34E-03 

0.01 

99.99 

TFORM 

4.42E-03 

X  4.42E-03 

0.01 

100.00 

Totals 

5;13E+01 

359504 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  'IN-LINE'  FACTOR  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 

(Factors  Greater  Than  1  Could  Indicate  Candidates  for  In-Lining) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time 

Percentage 

"In-Line" 

ARTCOMPE 

8.09E-01 

52400 

1.54E-05 

1.58 

28.84 

EVECTORE 

9.69E-01 

52400 

1.85E-05 

1.89 

24.08 

ARTCOMPZ 

1 . 18E+00 

34800 

3.38E-05 

2.29 

8.74 

EVECTORZ 

1 . 19E+00 

34800 

3.41E-05 

2.31 

8.68 

BETATVDE 

4.25E+00 

52400 

8.11E-05 

8.28 

5.49 

GSOLVE 

4.40E+00 

52400 

8.40E-05 

8.57 

5.30 

Factor 
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BETATVDZ 

4.69E+00 

34800  1.35E-04 

9.13 

2.19 

ESOLVE 

6.13E+00 

34800  1.76E-04 

11.94 

1.68 

BNDIN 

1.32E-02 

1000  1.32E-05 

0.03 

0.64 

BNDEX 

4.17E-02 

1000  4.17E-05 

0.08 

0.20 

BNDBLD 

1.97E-01 

1000  1.97E-04 

0.38 

0.04 

VISFLUX 

7.66E+00 

1000  7.66E-03 

14.92 

0.00 

ROEAVGE 

2.67E+00 

600  4.45E-03 

5.20 

0.00 

ROEAVGZ 

2.66E+00 

600  4.44E-03 

5.19 

0.00 

ATNSCl 

1.92E+00 

1  1.92E+00 

3.74 

0.00 

ALPHAZ 

1.61E+00 

600  2.68E-03 

3.13 

0.00 

EULFLUX 

1.59E+00 

1000  1.59E-03 

3.11 

0.01 

GCALCZ 

1.53E+00 

600  2.55E-03 

2.98 

0.00 

OUTPUT 

1.49E+00 

1  1.49E+00 

2.90 

0.00 

PSIZETA 

1.20E+00 

600  2.00E-03 

2.34 

0.00 

ALPHAE 

1.08E+00 

400  2.71E-03 

2.11 

0.00 

THSTEP 

i,08E+00 

200  5.38E-03 

2.09 

0.00 

GCALCE 

1.02E+00 

400  2.55E-03 

1.99 

0.00 

PSIETA 

8.05E-0i 

400  2.01E-03 

1.57 

0.00 

NORM 

5.01E-01 

80  6.26E-03 

0.98 

0.00 

EVALUEZ 

3.24E-01 

600  5.39E-04 

0.63 

0.01 

EVALUEE 

3.23E-01 

600  5.39E-04 

0.63 

0.01 

STORE 

7.56E-03 

20  3.78E-04 

0.01 

0.00 

INITIAL 

5.34E-03 

1  5.34E-03 

0.01 

0.00 

TFORM 

4.42E-03 

1  4.42E-03 

0.01 

0.00 

Totals 

5.13E+01 

359504 
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C.3  ATNSCl  (Anisotropic  and  Isotropic  Damping) 

The  data  processing  rate  for  the  variable  e  case  is  2.0457  x  10“®  seconds  per  grid 
point  per  time  level  for  the  CRAY  X-MP/216,  utiiling  a  133x60  grid.  FLOWTRA  A 
results  are  for  200  iterations  (400  time  levels). 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  TIME  USED  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

Ac  cum'/, 

BETATVDZ 

1.58E+01 

34800  4.54E-04 

24.22 

24.22 

VISFLUX 

7.63E+00 

1000  7.63E-03 

11.70 

35.91 

BETATVDE 

6.77E+00 

52400  1.29E-04 

10.37 

46.28 

FSOLVE 

6.13E+00 

34800  1.76E-04 

9.39 

55.67 

GSOLVE 

4.58E+00 

52400  8.74E-05 

7.02 

62.69 

ROEAVGE 

2.66E+00 

600  4.43E-03 

4.07 

66.77 

ROEAVGZ 

2.65E+00 

600  4.42E-03 

4.06 

70.83 

ATNSCl 

1.98E+00 

1  1.98E+00 

3.03 

73.86 

ALPHAZ 

1.60E+00 

600  2.66E-03 

2.45 

76.30 

EULFLUX 

1.58E+00 

1000  1.58E-03 

2.43 

73.73 

OUTPUT 

1.58E+00 

1  1.58E+00 

2.42 

81.15 

GCALCZ 

1.53E+00 

600  2.55E-03 

2.34 

83.49 

PSIZETA 

1.22E+00 

600  2.04E-03 

1.88 

85.37 

EVECTORZ 

1 . 18E+00 

34800  3.40E-05 

1.81 

87.18 

ARTCOMPZ 

1.16E+00 

34800  3.34E-05 

1.78 

88.96 

TMSTEP 

1.08E+00 

200  5.40E-03 

1.65 

90.61 

ALPHAE 

1.07E+00 

400  2.68E-03 

1.65 

92.26 

GCALCE 

1.02E+00 

400  2.55E-03 

1.56 

93.82 
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EVECTORE 

9.87E-01 

52400  1.88E-05 

1.51 

95.34 

PS I ETA 

8.16E-01 

400  2.04E-03 

1.25 

96.59 

ARTCOMPE 

8.01E-01 

52400  1.53E-05 

1.23 

97.81 

NORM 

5.16E-01 

84  6.15E-03 

0.79 

98.60 

EVALUEZ 

3.21E-01 

600  5.35E-04 

0.49 

99.10 

EVALUEE 

3.20E-01 

600  5.33E-04 

0.49 

99.59 

BNDBLD 

1.96E-01 

1000  1.96E-04 

0.30 

99.89 

BNDEX 

4.26E-02 

1000  4.26E-05 

0.07 

99.95 

BNDIN 

1.42E-02 

1000  1.42E-05 

0.02 

99.97 

STORE 

7.79E-03 

20  3.90E-04 

0.01 

99.99 

INITIAL 

5.37E-03 

1  5.37E-03 

0.01 

99.99 

TFORM 

4.35E-03 

1  4.35E-03 

0.01 

100.00 

Totals 

6.53E+01 

359508 

FLOWTRACE  RESULTS  OF 

ROUTINES 

SORTED  BY 

'IN-LINE'  FACTOR  (DESCENDING) 

(CPU  Times 

are  Shown  in 

Seconds) 

(Factors  Greater  Than  1  Could  Indicate  Candidates 

for  In-Lining) 

Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

"In-Line" 

ARTCOMPE 

8.01E-01 

52400  1.53E-05 

1.23 

29.11 

EVECTORE 

9.87E-01 

52400  1.88E-05 

1.51 

23.64 

ARTCOMPZ 

1.16E+00 

34800  3.34E-05 

1.78 

8.86 

EVECTORZ 

1 . 18E+00 

34800  3.40E-05 

1.81 

8.71 

GSOLVE 

4.58E+00 

52400  8.74E-05 

7.02 

5.09 

BETATVDE 

6.77E+00 

52400  1.29E-04 

10.37 

3.45 
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FSOLVE 

6.13E+00 

34800  1.76E-04 

9.39 

1.68 

BETATVDZ 

1.58E+01 

34800  4.54E-04 

24.22 

0.65 

BNDIN 

1.42E-02 

1000  1.42E-05 

0.02 

0.60 

BNDEX 

4.26E-02 

1000  4.26E-05 

0.07 

0.20 

BNDBLD 

1.96E-01 

1000  1.96E-04 

0.30 

0.04 

VISFLUX 

7.63E+00 

1000  7.63E-03 

11.70 

0.00 

ROEAVGE 

2.66E+00 

600  4.43E-03 

4.07 

0.00 

ROEAVGZ 

2.65E+00 

600  4.42E-03 

4.06 

0.00 

ATNSCl 

1.98E+00 

1  1.98E+00 

3.03 

0.00 

ALPHAZ 

1.60E+00 

600  2.66E-03 

2.45 

0.00 

EULFLUX 

1.58E+00 

1000  1.58E-03 

2.43 

0.01 

OUTPUT 

1.58E+00 

1  1.58E+00 

2.42 

0.00 

GCALCZ 

1.53E+00 

600  2.55E-03 

2.34 

0.00 

PSIZETA 

1.22E■^00 

600  2.04E-03 

1.88 

0.00 

TMSTEP 

1.08E+00 

200  5.40E-03 

1.65 

0.00 

ALPHAE 

1.07E+00 

400  2.68E-03 

1.65 

0.00 

GCALCE 

1.02E+00 

400  2.55E-03 

1.56 

0.00 

PS I ETA 

S.16E-01 

400  2.04E-03 

1.25 

0.00 

NORM 

5.16E-01 

84  6.15E-03 

0.79 

0.00 

EVALUEZ 

3.21E-01 

600  5.35E-04 

0.49 

0.01 

EVALUEE 

3.20E-01 

600  5.33E-04 

0.49 

0.01 

STORE 

7.79E-03 

20  3.90E-04 

0.01 

0.00 

INITIAL 

5.37E-03 

1  5.37E-03 

0.01 

0.00 

TFORM 

4.35E-03 

1  4.35E-03 

0.01 

0.00 

Totals 

6.53E+01 

359508 
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C.Ji  ATNSCS  (No  Jacobian  Update  Between  Operators) 


The  data  processing  rate  when  the  viscous  Jacobians  are  updated  onl}-  after 
a  complete  sequence  of  operator  sweeps  is  7.6128  x  10“^  seconds  per  grid  point  per 
time  level.  This  is  for  the  CRAY  X-MP/216,  utiiling  a  133  x  60  grid.  FLOVVTR.ACE 
results  are  for  200  iterations  (400  time  levels). 


FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  TIME  USED  (DESCENDING) 
(CPU  Times  are  Shown  in  Seconds) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

Accum'/, 

JACOBIAN 

1.61E+02 

200  8.05E-01 

66.26 

66.26 

PSIZETA 

2.01E+01 

600  3.35E-02 

8.26 

74.53 

PSIETA 

1.34E+01 

400  3.36E-02 

5.53 

80.05 

VISFLUX 

7.52E+00 

1000  7.52E-03 

3.10 

83.15 

FSOLVE 

6.03E+00 

34800  1.73E-04 

2.48 

85.63 

BETATVDZ 

4.69E+00 

34800  1.35E-04 

1.93 

87.56 

GSOLVE 

4.34E+00 

52400  8.28E-05 

1.79 

89.54 

BETATVDE 

3.89E+00 

52400  7.43E-05 

1.60 

90.95 

ROEAVGE 

2.66E+00 

600  4.43E-03 

1.09 

92.04 

ROEAVGZ 

2.65E+00 

600  4.41E-03 

1.09 

93.13 

ATNSC2 

1.92E+00 

1  i.92E+00 

0.79 

93.92 

ALPHAZ 

1.58E+00 

600  2.63E-03 

0.65 

94.57 

EULFLUX 

1.57E+00 

1000  1.57E-03 

0.65 

95.21 

OUTPUT 

1 . 54E+00 

1  1.54E+00 

0.63 

95.85 

GCALCZ 

1.52E+00 

600  2.54E-03 

0.63 

96.47 

EVECT.ORZ 

1.16E+00 

34800  3.34E-05 

0.48 

96.95 

ARTCOMPZ 

1.14E+00 

34800  3.27E-05 

0.47 

97.42 
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TMSTEP 

1.07E+00 

200  5.36E-03 

0.44 

97.86 

ALPHAS 

1.06E+00 

400  2.66E-03 

0.44 

98.30 

GCALCE 

1.02E+00 

400  2.54E-03 

0.42 

98.72 

EVECTORE 

9.42E-01 

52400  1.80E-05 

0.39 

99.10 

ARTCOMPE 

7.77E-01 

52400  1.48E-05 

0.32 

99.42 

NORM 

5.00E-01 

80  6.26E-03 

0.21 

99.63 

EVALUEZ 

3.14E-01 

600  5.23E-04 

0.13 

99.76 

EVALUEE 

3.11E-01 

600  5.19E-04 

0.13 

99.89 

BNDBLD 

2.03E-01 

1000  2.03E-04 

0.08 

99.97 

BNDEX 

4.14E-02 

1000  4.14E-05 

0.02 

99.99 

BNDIN 

1.31E-02 

1000  1.31E-05 

0.01 

99.99 

STORE 

7.45E-03 

20  3.73E-04 

0.00 

100.00 

INITIAL 

5.40E-03 

1  5.40E-03 

0.00 

100.00 

TFORti 

4.22E-03 

1  4.22E-03 

0.00 

100.00 

Totals 


2.43E+02  359704 


FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  'IN-LINE’  FACTOR  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 

(Factors  Greater  Than  1  Could  Indicate  Candidates  for  In-Lining) 


Routine  Name 

Tot  Time  # 

Calls  Avg  Time  Percentage 

"In-Line"  Factor 

ARTCOMPE 

7.77E-01 

52400  1.48E-05 

0.32 

30.03 

EVECTORE 

9.42E-01 

52400  1.80E-05 

0.39 

24.76  ■ 

ARTCOMPZ 

1.14E+00 

34800  3.27E-05 

0.47 

9.05 

EVECTORZ 

1.16E+00 

34800  3.34E-05 

0.48 

8.86 

BETATVDE 

3.89E+00 

52400  7.43E-05 

1.60 

5.99 

GSOLVE 

4.34E+00 

52400  8.28E-05 

1.79 

5.38 

BETATVDZ 

4.69E+00 

34800  1.35E-04 

1.93 

2.19 

FSOLVE 

6.03E+00 

34800  1.73E-04 

2.48 

1.71 

BNDIN 

1.31E-02 

1000  1.31E-05 

0.01 

0.65 

BNDEX 

4.14E-02 

1000  4.14E-05 

0.02 

0.21 

BNDBLD 

2.03E-01 

1000  2.03E-04 

0.08 

0.04 

JACOBIAN 

1 .61E+02 

200  8.05E-01 

66.26 

0.00 

PSIZETA 

2.01E+01 

600  3.35E-02 

8.26 

0.00 

PS I ETA 

1.34E+01 

400  3.36E-02 

5.53 

0.00 

VISFLUX 

7.52E+00 

1000  7.52E-03 

3.10 

0.00 

ROEAVGE 

2.66E+00 

600  4.43E-03 

1.09 

0.00 

ROEAVGZ 

2.65E+00 

600  4.41E-03 

1.09 

0.00 

ATNSC2 

1.92E+00 

1  1.92E+00 

0.79 

0.00 

ALPHAZ 

1.58E+00 

600  2.63E-03 

0.65 

0.00 

EULFLUX 

1.57E+00 

1000  1.57E-03 

0.65 

0.01 

C-139 


oyiPUT 

1.54E+00 

1  1.54E+00 

0.63 

0.00 

GCALCZ 

1.52E+00 

600  2.54E-03 

0.63 

0.00 

TMSTEP 

1.07E+00 

200  5.36E-03 

0.44 

0.00 

ALPHAE 

1.06E+00 

400  2.66E-03 

0.44 

0.00 

GCALCE 

1.02E+00 

400  2.54E-03 

0.42 

0.00 

NORM 

5.00E-01 

80  6.26E-03 

0.21 

0.00 

EVALUEZ 

3.14E-01 

600  5.23E-04 

0.13 

0.01 

EVALUEE 

3.11E-01 

600  5.19E-04 

0.13 

0.01 

STORE 

7.45E-03 

20  3.73E-04 

0.00 

0.00 

INITIAL 

5.40E-03 

1  5.40E-03 

0.00 

0.00 

TFORM 


4.22E-03 


1  4.22E-03 


0.00 


0.00 


G.o  ATNSC2  (Jacobian  Update  After  Each  Operator  Sweep) 


The  data  processing  rate  when  the  viscous  Jacobians  are  updated  after  each 
operator  sweep  is  1.9988  x  lO”"*  seconds  per  grid  point  per  time  level.  This  is  for 
the  CRAY  X-MP/216,  utiiling  a  133  x  60  grid.  FLOVVTRACE  results  are  for  200 
iterations  (400  time  levels). 

FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  TIME  USED  (DESCENDING) 

(CPU  Times  are  Shown  in  Seconds) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

Accum*/, 

PSIZETA 

2.41E+02 

600  4.01E-01 

37.77 

37.77 

JACOBIAN 

1.86E+02 

1000  1.86E-01 

29.15 

66.93 

PSIETA 

1.62E+02 

400  4.04E-01 

25.37 

92.30 

VISFLUX 

7.54E+00 

1000  7.54E-03 

1.18 

93.48 

FSOLVE 

6.09E+00 

34800  1.75E-04 

0.95 

94.44 

BETATVDZ 

4.69E+00 

34800  1.35E-04 

0.74 

95.17 

GSOLVE 

4.56E+00 

52400  8.69E-06 

0.71 

95.89 

BETATVDE 

4.16E+00 

52400  7.93E-05 

0.65 

96.54 

ROEAVGE 

2.65E+00 

600  4.42E-03 

0.42 

96.95 

ROEAVGZ 

2.65E+00 

600  4.41E-03 

0.41 

97.37 

ATNSC2 

1.93E+00 

1  1.93E+00 

0.30 

97.67 

ALPHAZ 

1.58E+00 

600  2.63E-03 

0.25 

97.92 

EULFLUX 

1.57E+00 

1000  1.57E-03 

0.25 

98.17 

OUTPUT 

1.55E+00 

1  1.55E+00 

0.24 

98.41 

GCALCZ 

1.53E+00 

600  2.55E-03 

0.24 

98.65 

EVECTORZ 

1.16E+00 

34800  3.33E-05 

0.18 

98.83 

ARTCOMPZ 

1 . 14E+00 

34800  3.29E-05 

0.18 

99.01 

TMSTEP 

1.06E+00 

200 

5.32E-03 

0.17 

99.18 

ALPHAE 

1.06E+00 

400 

2.66E-03 

0.17 

99.34 

GCALCE 

1.02E+00 

400 

2.54E-03 

0.16 

99.50 

EVECTORE 

9.79E-01 

52400 

1.87E-05 

0.15 

99.66 

ARTCOMPE 

7.85E-01 

52400 

1.50E-05 

0.12 

99.78 

NORM 

5.03E-01 

80 

6.29E-03 

0.08 

99.86 

EVALUEZ 

3.14E-01 

600 

5.24E-04 

0.05 

99.91 

EVALUEE 

3.13E-01 

600 

5.22E-04 

0.05 

99.96 

BNDBLD 

2.03E-01 

1000 

2.03E-04 

0.03 

99.99 

BNDEX 

4.24E-02 

1000 

4.24E-05 

0.01 

100.00 

BNDIN 

1.42E-02 

1000 

1.42E-05 

100.00 

STORE 

7.44E-03 

20 

3.72E-04 

100.00 

INITIAL 

5.35E-03 

1 

5.35E-03 

100.00 

TFORM 

4.32E-03 

1 

4.32E-03 

100.00 

Totals 


6.38E+02  360504 


FLOWTRACE  RESULTS  OF  ROUTINES 
SORTED  BY  'IN-LINE'  FACTOR  (DESCENDING) 

(CPU  Times<are  Shown  in  Seconds) 

(Factors  Greater  Than  1  Could  Indicate  Candidates  for  In-Lining) 


Routine  Name 

Tot  Time 

#  Calls  Avg  Time  Percentage 

"In-Line" 

ARTCOMPE 

7.85E-01 

52400  1.50E-05 

0.12 

29.73 

EVECTORE 

9.79E-01 

52400  1.87E-05 

0.15 

23.82 

ARTCOMPZ 

1.14E+00 

34800  3.29E-05 

0.18 

8.99 

EVECTORZ 

1 . 16E+00 

34800  3.33E-05 

0.18 

8.87 

BETATVDE 

4.16E+00 

52400  7.93E-05 

0.65 

5.61 

GSOLVE 

4.56E+00 

52400  8.69E-05 

0.71 

5.12 

BETAtVDZ 

4.69E+00 

34800  1.35E-04 

0.74 

2.19 

FSOLVE 

6.09E+00 

34800  1.75E-04 

0.95 

1.69 

BNDIN 

1.42E-02 

1000  1.42E-05 

0.00 

0.60 

BNDEX 

4.24E-02 

1000  4.24E-05 

0.01 

0.20 

BNDBLD 

2.03E-01 

1000  2.03E-04 

0.03 

0.04 

PSIZETA 

2.41E+02 

600  4.01E-01 

37.77 

0.00 

JACOBIAN 

1.86E+02 

1000  1.86E-01 

29.15 

0.00 

PS I ETA 

1.62E+02 

400  4.04E-01 

25.37 

0.00 

VISFLUX 

7 . 54E+00 

1000  7.54E-03 

1.18 

0.00 

ROEAVGE 

2.65E+00 

600  4.42E-03 

0.42 

0.00 

ROEAVGZ 

2.65E+00 

600  4.41E-03 

0.41 

0.00 

ATNSC2 

1.93E+00 

1  1.93E+00 

0.30 

0.00 

ALPHAZ 

1 . 58E+00 

600  2.63E-03 

0.25 

0.00 

EULFLUX 

1 . 57E+00 

1000  1.57E-03 

0.25 

0.01 

Factor 


C-I'i:] 


OUTPUT 

•  1.55E+00 

1  1.55E+00 

0.24 

0.00 

GCALCZ 

1 . 53E+00 

600  2.55E-03 

0.24 

0.00 

TMSTEP 

1.06E+00 

200  5.32E-03 

0.17 

0.00 

ALPHAE 

1.06E+00 

400  2.66E-03 

0.17 

0.00 

GCALCE 

1.02E+00 

400  2.54E-03 

0.16 

0.00 

NORM 

5.03E-01 

80  6.29E-03 

0.08 

0.00 

Evaluez 

3.14E-01 

600  5.24E-04 

0.05 

0.01 

EVALUEE 

3.13E-01 

600  5.22E-04 

0.05 

0.01 

STORE 

7.44E-03 

20  3.72E-04 

0.00 

0.00  ■ 

INITIAL 

5.35E-03 

1  5.35E-03 

0.00 

0.00 

TFORM 

4.32E-03 

1  4.32E-03 

0.00 

0.00 

totals 

6.38E+02 

360504 
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